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ABSTRACT 

We determine the neutral kaon mixing matrix element Bk in the continuum limit with 2+1 flavors 
of domain wall fermions, using the Iwasaki gauge action at two different lattice spacings. These 
lattice fermions have near exact chiral symmetry and therefore avoid artificial lattice operator 
mixing. 

We introduce a significant improvement to the conventional NPR method in which the bare ma- 
trix elements are renormalized non-perturbatively in the RI-MOM scheme and are then converted 
into the MS scheme using continuum perturbation theory. In addition to RI-MOM, we intro- 
duce and implement four non-exceptional intermediate momentum schemes that suppress infrared 
non-perturbative uncertainties in the renormalization procedure. We compute the conversion fac- 
tors relating the matrix elements in this family of RI-SMOM schemes and MS at one-loop order. 
Comparison of the results obtained using these different intermediate schemes allows for a more 
reliable estimate of the unknown higher-order contributions and hence for a correspondingly more 
robust estimate of the systematic error. We also apply a recently proposed approach in which 
twisted boundary conditions are used to control the Symanzik expansion for off-shell vertex func- 
tions leading to a better control of the renormalization in the continuum limit. 
We control chiral extrapolation errors by considering both the NLO SU(2) chiral effective theory, 
and an analytic mass expansion. We obtain 5^^(3GeV) = 0.529(5)stat(15);f;(2)Fv(l 1)npr- This 
corresponds to B\^^ = 0.749(7)stat(21);^(3)Fv(15) NPR. Adding all sources of error in quadrature 
we obtain = 0.749(27) 

combined J with an Overall combined error of 3.6%. 
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I. INTRODUCTION 

The indirect CP violation parameter of the neutral kaon system 

_ A{KL^{mt)i=Q) 
A{Ks ^ {7t7t)i=oy 

n 

was measured first at BNL in a Nobel Prize winning experiment [1], and is now experimentally 
measured as \£k\ = (2.228 ± 0.01 1) 10^^ ||2|]. Since CP is not an exact symmetry of the weak 
interations, the eigenstates Ki and ^5 of the mass matrix of neutral kaon system are not eigenstates 
of CP We characterise the state mixing via 

Ks = pK^-qK° and KL = pK^ + qK^ (2) 

where + = I, and ^ = j^. 

£k receives its dominant contribution from "indirect" CP violation via state-mixing, mediated by 
the imaginary part of the AS = 2 box graph. Before £k can be used to constrain the unitarity 
triangle and to provide information on CKM matrix elements, we must therefore determine the 
QCD hadronic matrix element of the effective weak A5 = 2 four quark operator 



(^°|(^vv+aa|^), 
where 

^vv+AA = {sy^d) {sy^d) + [sy^y^d) [sy^y^d) . (3) 
It is conventional to define the bag parameter Bk from this matrix element as 

where Mk and fK are the mass and leptonic decay constant of the kaon. The kaon bag parameter 
is thus of fundamental importance in studies of CP violation, and as the hadronic matrix element 
is non-perturbative, lattice QCD is the only known framework for its determination from first 
principles. 

Since the operator ^vv+aa depends on the renormalization scheme and scale used in its definition, 
Bk also has the same scheme and scale dependence. Therefore, for phenomenological use, it is 
convenient to introduce the renormalization-group-invariant counterpart of Bk, 

BK = (OA\n,nf)B^{^,nf), 
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where the Wilson coefficient, CO^^ {ji.nf), for the various schemes A used in this paper are given 
in Equations (|66l) through (TTOl) . and we use the numerical values for the 2+1 flavour theory in our 
conversion. 

We have recently calculated Bk in dynamical 2+1 flavored simulations [|3L |4|]with a total error 
of about 5.5%. It was observed by Buras and Guadagnoli [|5|], that our result f3\ was sufficiently 
accurate that additional care needs to be taken in relating it to the measured value of Ek- Previously 
ignored subdominant effects of direct CP violation arising from the A5 = 1 Hamiltonian amount 
to a few percent and must now be incorporated. 

The short distance contribution Ca" 7] differs from £k, predominantly due to direct CP violation 



£k = £k + i 



.ImAo 



ReAo' 



Here Aq is the — >■ nn amplitude for the isospin final state defined via 



A{K^ ^7t7t{I))=Aiexpidi and A ^ 7t7t{I)) = A*i expi5i 



(5) 



(6) 



and 5/ is the ktt phase shift in the / = or / = 2 final state. 

Reliable calculation of Aq amplitudes remains a challenging project to which our collaboration is 
devoting a considerable effort [|8l-ll3n. Using the measured value = (1.65±0.26) x 10^ |12|], 
assuming the Standard Model is correct and making plausible assumptions in estimating the some- 
what less difficult ratio ^^i, the subdominant contribution to £k can be effectively incorporated 
into a correction factor fCg^ [|5D: 



(7) 



where Xx = V^dV*^ contain the entries of the CKM matrix Vxy, rji are perturbative QCD corrections 



n 



151] the cor- 



[|14|] and the Sq are Inami-Lim functions of mass ratios Xq = —r. In References p, 
rection factor was estimated to be ^ 0.94 ± 0.02, and here the fractional error on this small 
correction is large (0.02 in a correction of size 0.06) and model dependent. 
The correction factor also includes an estimate of long distance contributions corresponding to 
two insertions of the A5 = 1 Hamiltonian, with two pions propagating long distances between 
them [|l5ll . The results of our present work are sufficiently precise that it has become necessary to 
determine as many contributions as possible using lattice gauge methods; efforts in RBC-UKQCD 



are underway in this direction 
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In this paper we improve on our earlier calculations in three major ways. First of all, we 

simulate at a second value of the lattice spacing which allows us to perform a continuum ex- 
trapolation. Secondly, we refine our approach to non-perturbative renormalization to implement 
intermediate schemes defined with no exceptional momentum channels and thereby reduce the in- 
frared non-perturbative uncertainties. Finally, we also use twisted boundary conditions to remove 
the requirement to use the Fourier modes of our lattice for our renormalization of off- shell ampli- 
tudes: this gives complete freedom of choice of the momentum at each lattice spacing and enables 
a more reliable continuum extrapolation of the renormalized operator. 

Our final result for Bk from the present analysis is obtained using an off- shell momentum scheme 
renormalization. When converted to MS with = fi^ = (3GeV)~ it is: 



5^S(3GeV)=0.529(5)stat(15);^(2)Fv(ll) NPR- 



(8) 



The 3 GeV scale for our result is made accessible by our improved renormalization techniques, and 
enables us to reduce perturbative error compared to a 2 GeV renormalization scale. For comparison 
to other results we also quote the standard operator normalization: 



= o.749(7)stat(21);t(3)Fv(15)NPR. 



(9) 



The full analysis of systematic errors presented in this paper augments and finalizes an earlier 



conference presentation lllSh . The result Equation ([8]) represents around a factor of four reduction 
in the error during the last five years or so. 

The structure of the remainder of this paper is as follows. In the next section we discuss the details 
of our simulations and present the measured values of the bare matrix elements. In Section Hill 
we discuss the definition of several new momentum renormalization schemes and perform the 
non-perturbative renormalization of the bare lattice operator i^?'vv+aa into these schemes. In this 
section we also perform the one-loop perturbative matching from the momentum schemes into MS. 
Having obtained the matrix elements at the values of the quark masses and lattice spacing at which 
we perform our simulations, we present the simultaneous chiral and continuum extrapolations of 
the renormalized matrix elements in Section |IVl We will discuss the phenomeno logical context of 
our results in the concluding Section |VI] of this paper. 
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Lattice 



1 (32^ X 64) 



2 (24^ X 64) 



rrih mi 



traj.(# meas.) 



0.03 0.004 260-3250 (300) 
0.03 0.006 500-3610 (312) 
0.03 0.008 260-2770 (252) 



0.04 0.005 900-8940 (202) 
0.04 0.01 1460-8540(178) 

TABLE I: Ensemble details. Here traj. refers to the Monte Carlo trajectories used in our measurements. The 
bracketed # meas. refers to the number of measurements, separated by 20 MD time units (10 trajectories) 
for the 1 ensembles, and 40 molecular- dynamics time units (40 trajectories) for the 2 ensembles. To 
reduce the effects of auto-correlations we block-average our data over 80 MD time units and use blocked 
measurements for the purposes of statistical analysis. 

n. SIMULATION PARAMETERS AND MATRIX ELEMENTS 



Details of our ensembles are given in references 



Iwasaki gauge action 12011 with 2-1-1 flavors of dynamical domain wall fermions 



1911 . and are summarised in Table Jl We use the 

. This action 



was chosen to balance topology change against chirality after a careful study |122 



21 



2411 recognising 



a general problem that topological tunneling will vanish towards the continuum limit in any local 
update due to the gauge field potential barrier 11221 



2511 ■ These lattice fermions have near exact 



chiral symmetry and avoid artificial lattice operator mixing, while retaining acceptable topology 
change in our region of simulation. 

We have two lattices of similar physical volume at two lattice spacings: 

(i) Our finer lattice has 32^ x 64 x 16 points and a coupling /3 = 2.25, which our analysis sug- 
gests corresponds to an inverse lattice spacing a^^ = 2.28(3) GeV. We refer to the ensembles 
with /3 — 2.25 as the 1 ensemble set. 

(ii) Our coarser lattice has 24^ x 64 x 16 points and a coupling /3 = 2.13, corresponding to 
a^^ = 1.73(3) GeV. The ensembles with /3 = 2.13 are labeled as the 2 ensemble set. 

For each ensemble set we use a number of valence masses to increase the amount of information 
in the light mass regime. We use our standard notation for quark masses, m/ and m^, represent 
respectively the lighter and heavier of the two sea-quark masses (the sea consists of two quarks 
with mass m/ and one with mass m/,). For the valence masses we use subscripts from the end 
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Lattice 


nih 


{m/} 




1 32^ X 64 


0.03 


0.004,0.006,0.008 


0.002, 0.004, 0.006, 0.008, 0.025, 0.03 


2 24^ X 64 


0.04 


0.005,0.01 


0.001,0.005,0.01,0.03,0.02,0.04 



TABLE II: Details of partially quenched valence masses {m,,} on each ensemble. Meson correlation func- 
tions were computed for all possible pairings of valence masses. 

of the alphabet my, and niy as appropriate. /, are masses in the DWF action used in the 
simulation whereas the valence masses appear in the corresponding partially quenched action. 
Because of the finite extent of the fifth dimension, small residual mass effects are present and 
the multiplicatively renormalizable bare quark masses are defined as m/,/,,v,x,.y = ^i,h,v,x,y-\-^Ks-, 
where mj-es is the residual mass. The values of the valence quark masses used in our measurements 
are summarised in Table HIl As in Reference we will restrict our analysis, which relies on 
SU(2) chiral perturbation theory, to light-quark masses corresponding to pions lighter than about 
420 MeV. 

We use two approaches to calculate the matrix element (^''l ^vv+aaI-^'^)- Both combine periodic 
and anti-periodic boundary conditions in the time direction to eliminate the leading, unwanted 
around-the-world propagation of the meson states that arise with a finite lattice in the time direc- 
tion. In both cases we use gauge-fixed wall sources to create a state and annihilate a state, 
and form a ratio 

{K\t,)Wvv+AA{t)\K^{t2)) 



B 



lat 
K 



(10) 



l{KOin)\Ao{t)){Ao{t)m2)) 
For convenience we use the local axial current interpolating operators in the denominator, and this 
ratio must be multiplied by a renormalization constant 



'-^VV+AA 



(11) 



to obtain physically normalized matrix elements. 

On our 1 ensembles we used a single source at ? = and used the {P+A) combination for the 
forward propagating K meson. This has the effect of creating {P +A) x (P + A) = PP + AA + 
PA-\-AP combinations in meson propagators, and the meson state has periodicity 2Lj, where 
Lj- = 64 is the temporal extent of the lattice. Similarly the {P —A) combination is taken for the 
backward propagating K meson. These Fermion boundary conditions are implemented on gauge 
links crossing the toroidal wrapping plane between t = and r = — 1 . On each successive gauge 
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0.134 



0.132 



0.13 



' 0.128 ■ 



0.126 



0.124 



0.122 



12 16 20 24 28 32 36 40 44 48 52 56 
t 



FIG. 1: Effective mass plateau of the lightest unitary simulated pion (m/, = 0.03, = ruy = mi = 0.004) on 
the 1 ensembles. Here the plateau is obtained from the wall-local PP correlator, but the fit displayed is to 
all pseudoscalar correlators. 



configuration we selected a different time ?src at which to insert the kaon sources. For simplicity 
this was implemented by translating the gauge configuration and redefining /^src to be zero. The 
boundary condition described above is then applied. 

The above approach requires half the number of propagator inversions on each configuration (and 
enables us to sample more frequently at fixed cost) compared to that taken on the 2 ensembles. 
On our 2 ensembles we used a source at ? = 5 and a source dXt = 59 requiring seperate inversions 
for each source. For each propagator entering a meson, we took the average of periodic (P) and 
anti-periodic (A) solutions. 

The A5 = 2 four-quark operator i^vv+aa is inserted on all times between the kaon creation and 
anti-kaon annihilation operators. The locations of the kaon, anti-kaon and operator all receive 1? 
volume averages, giving a low variance estimate of the correlation function. 
The quality of the data can be gauged from Figures [T] through [6l displaying the lightest simulated 
pion, heaviest eta and a typical kaon matrix element fit to B^^ for each of the two lattice spacings. 
More examples can be found in ref lll9n . Tables Hill and ITVl display the fitted values for the matrix 



element on each lattice. The fitted meson masses are as in reference 
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FIG. 2: Effective mass plateau of the lightest unitary simulated pion (m/, = 0.04, nix = = "J/ = 0.005) on 
the 2 ensembles. Here the plateau is obtained from the wall-local PP correlator, but the fit displayed is to 
all pseudoscalar correlators. 



0.33 



0.328 - 



0.326 



0.324 



0.322 - 



0.32 



0.318 



0.316 




FIG. 3: Effective mass plateau of the heaviest simulated eta {nix = = "ift = 0.03, m/ = 0.008) on the 
1 ensembles. Here the plateau is obtained from the wall-local PP correlator, but the fit displayed is to all 
pseudoscalar correlators. 



A. Reweighting 



As explained above, at each lattice spacing we have performed the simulations using a number of 
light-quark masses but only a single sea strange-quark mass. As we can only determine the phys- 
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FIG. 4: Effective mass plateau of the heaviest simulated eta (nix = "^y = "^/i = 0.04, nii = 0.01) on the 
2 ensembles. Here the plateau is obtained from the waU-local PP correlator, but the fit displayed is to all 
pseudoscalar correlators. 
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FIG. 5: A typical B^' matrix element correlator (niy = = 0.03, = mi = 0.004) on the 1 ensembles. 



ical strange quark mass after the analysis is complete, our imperfect pre-simulation estimate 
of Ms has been a source of error in previous calculations, where we could only adjust the valence 
strange quark mass or use SU(3) chiral perturbation theory to estimate the effects of varying the 
unitary strange quark mass. We do not expect significant effects from small adjustments of the sea 
strange-quark mass and reweighting gives us a tool to demonstrate this without doubling the cost 



of the simulation. For more discussion we refer to our papers [|4, 



m. 
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FIG. 6: A typical B^' matrix element correlator {niy = m/, = 0.04, = mi = 0.005) on the 2 ensembles. 



0.62 



0.61 



0.59 



■ 0.58 




0.57 



0.56 



0.55 - 



0.54 



FIG. 7: An overlay of a typical 5^' matrix element (niy = 0.03, = mj = 0.004) on the 1 ensembles at 
two values of the sea strange quark mass: = 0.03 (red) and m/, = 0.027 (blue). The latter is at our closest 
reweight to the physical strange mass. 

Figure |7] shows an overlay of a typical kaon matrix element correlator at the simulated sea 
strange-quark mass and the physical value. Figure [8] shows the dependence of the fitted value of 
the matrix element of ^vv+aa on the sea strange-quark mass; the dependence is very small and 
barely statistically significant. 

For both ensemble sets, we compute the propagators at two valence strange-quark masses: my ~ 
0.03 and 0.025 for the 1 ensembles and lUy = 0.04 and 0.03 for the 2 ensembles. When 



12 



My 


B^(m/ = 0.004) B 


^{mi = 0.006) B^{mi = 0.008) 


0.03 0.03 


0.6289(12) 


0.6305(12) 


0.6295(12) 


0.025 0.03 


0.6199(12) 


0.6214(12) 


0.6207(12) 


0.008 0.03 


0.5862(17) 


0.5878(17) 


0.5878(19) 


0.006 0.03 


0.5823(19) 


0.5838(21) 


0.5838(22) 


0.004 0.03 


0.5787(24) 


0.5801(27) 


0.5798(28) 


0.002 0.03 


0.5767(46) 


0.5772(43) 


0.5781(50) 


0.025 0.025 


0.6100(13) 


0.6116(13) 


0.6110(13) 


0.008 0.025 


0.5725(16) 


0.5745(17) 


0.5741(18) 


0.006 0.025 


0.5679(17) 


0.5701(20) 


0.5694(21) 


0.004 0.025 


0.5634(21) 


0.5659(24) 


0.5649(26) 


0.002 0.025 


0.5601(39) 


0.5629(37) 


0.5630(43) 


0.008 0.008 


0.5135(18) 


0.5178(19) 


0.5141(20) 


0.006 0.008 


0.5047(19) 


0.5096(20) 


0.5056(22) 


0.004 0.008 


0.4951(21) 


0.5013(23) 


0.4969(25) 


0.002 0.008 


0.4852(28) 


0.4939(32) 


0.4901(34) 


0.006 0.006 


0.4949(20) 


0.5004(22) 


0.4961(24) 


0.004 0.006 


0.4842(23) 


0.4908(25) 


0.4864(27) 


0.002 0.006 


0.4727(29) 


0.4813(34) 


0.4781(35) 


0.004 0.004 


0.4721(26) 


0.4791(29) 


0.4753(31) 


0.002 0.004 


0.4584(32) 


0.4663(37) 


0.4647(39) 


0.002 0.002 


0.4408(39) 


0.4473(44) 


0.4500(48) 



TABLE III: Fitted B^^ matrix element values on the 1 ensembles. For heavy-light matrix elements, my is 
the heavy quark mass. We chose a fit range of f = 12 — 52. 

computing kaonic quantities we reweight the sea strange mass irih to both valence strange-quark 
masses niy such that m/, = my in our observables. For each lattice and at each value of m/ we 
therefore have results with two strange quark masses with m/, = niy, one at the strange-quark mass 
at which we perform the simulation and the second obtained by reweighting. This enables us to 
interpolate linearly in the unitary strange quark mass to the physical point. In Tables IVl and IVTl we 
give the values for the heavy-light Bxy matrix element on each ensemble; it is to these data that we 
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My 


Bxy{m, = 0.005) B 


^(m;=0.01) 


0.04 


0.04 


0.6565(12) 


0.6562(12) 


0.03 


0.04 


0.6435(14) 


0.6430(13) 


0.02 


0.04 


0.6298(16) 


0.6291(14) 


0.01 


0.04 


0.6154(20) 


0.6145(17) 


0.005 


0.04 


0.6081(26) 


0.6078(24) 


0.001 


0.04 


0.6017(48) 


0.6072(53) 


0.03 


0.03 


0.6286(14) 


0.6280(13) 


0.02 


0.03 


0.6124(16) 


0.6117(14) 


0.01 


0.03 


0.5949(19) 


0.5943(16) 


0.005 


0.03 


0.5860(23) 


0.5860(20) 


0.001 


0.03 


0.5787(40) 


0.5835(40) 


0.02 


0.02 


0.5929(17) 


0.5924(15) 


0.01 


0.02 


0.5712(19) 


0.5711(16) 


0.005 


0.02 


0.5598(23) 


0.5603(19) 


0.001 


0.02 


0.5505(36) 


0.5547(31) 


0.01 


0.01 


0.5431(22) 


0.5439(18) 


0.005 


0.01 


0.5272(26) 


0.5284(21) 


0.001 


0.01 


0.5134(37) 


0.5164(29) 


0.005 0.005 


0.5075(31) 


0.5085(24) 


0.001 0.005 


0.4893(42) 


0.4903(31) 


0.001 0.001 


0.4652(55) 


0.4631(40) 



TABLE IV: Fitted matrix element values on the 2 ensembles. For heavy-light matrix elements, niy is 
the heavy quai^k mass. We chose a fit range of f = 12 — 52. 

perform our simultaneous chiral fits in Section ITVl 

III. NON-PERTURBATIVE RENORMALISATION 

In this section we discuss the renormalization of the A5 = 2 operator ^vv+aa? whose matrix 
elements we are computing. We start by performing non-perturbative renormalization, calculat- 
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FIG. 8: The m/, dependence of a typical B^' matrix element (my = 0.03, = mi = 0.004) on the 1 ensem- 
bles. 



nix 


Bxh{mi = 0.004) B.v/,(m; = 0.006) Bv/,(m, = 0.008) 


0.008 


0.5802(27) 


0.5807(29) 


0.5829(26) 


0.006 


0.5758(29) 


0.5764(32) 


0.5789(29) 


0.004 


0.5715(33) 


0.5721(38) 


0.5752(36) 


0.002 


0.5679(49) 


0.5680(52) 


0.5742(59) 



TABLE V: Heavy-light B^' matrix element values on the 1 ensembles at the physical ruh = Q.OTTiiJ), 
nih + m^efi = 0.0278(7) obtained from the NLO PQChPT combined fits of Section |Vl These values are 
obtained by first reweighting to m/, = niy then linearly interpolating in the unitary strange mass. 



rrix 


Bxh{mi =0.005) B 


,/,(m/ = 0.01) 


0.02 


0.6191(32) 


0.6190(27) 


0.01 


0.6035(35) 


0.6029(31) 


0.005 


0.5959(38) 


0.5949(37) 


0.001 


0.5892(64) 


0.5904(65) 



TABLE VI: Heavy-light B^ matrix element values on the 2 ensembles at the physical nt/, = 0.035(1), 
ruk + ?Mres = 0.038(1) obtained from the NLO PQChPT combined fits of Section |Vl These values ai^e 
obtained by first reweighting to m/, = ruy then linearly interpolating in the unitary strange mass. 
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ing numerically the renormalization factor which relates the bare lattice operator corresponding 
to our choice of the discrete QCD action to that defined in some intermediate renormalization 
scheme. For this to be feasible, of course, it is necessary that the intermediate scheme can be 
implemented numerically and we use several momentum subtraction schemes which are general- 
izations of the original Rl-MOM scheme i26|] . In phenomenological applications, our results for 
the matrix element (/^^|i^?'vv+aa|^) have to be combined with the Wilson coefficient function 
which is calculated in perturbation theory, most frequently using renormalization schemes based 
on dimensional regularization, such as the NDR scheme. It is therefore necessary to combine the 
coefficient function and the operator matrix element in the same scheme. Below we present the 
matching factors which relate the operator renormalized in our intermediate schemes to the corre- 
sponding operator in the NDR scheme. Since dimensional regularization cannot be implemented 
in lattice simulations, this (continuum) matching is performed in perturbation theory (at one-loop 
order) and is of course independent of the lattice calculations. The procedure described above can 
be summarised as follows: 

NPR 

Bare Lattice Operator — )■ Renormalized Operator in Momentum Subtraction Scheme 



Perturbation Theory , ^ . tttt h-tt^t^ n i 

Renormalized Operator in MS-NDR Scheme. 

The momentum subtraction schemes which we use require the evaluation of the Green functions 
for the transition d{pi)s{p2) — )■ d{p3)s{p4) with a suitable choice of the momenta p,. In the past. 



see in particular Reference Is], the results were presented using the Rl-MOM kinematic configu- 
ration in which pi = —p2 = P3 = —p4 1271] . Whilst this is correct asymptotically, i.e. when the pj 



3 



that 



are sufficiently large for each choice of the quark masses, it was argued in References [|28- 
performing the renormalization using Green functions with no exceptional channels, i.e. with no 
channels in which the square of the momentum is small, suppresses the non-asymptotic chiral 
symmetry breaking effects more effectively. In addition to the theoretical arguments, numerical 
evidence was presented demonstrating the suppression of terms which violated the chiral Ward- 
Takahashi identities, such as the equality of the renormalization constants of the vector and axial 
currents and of the scalar and pseudoscalar densities. Although the effects are small, typically of 
the order of a few percent, lattice calculations are becoming sufficiently precise that the reduction 
of this systematic error is necessary. 

For Bk, the Rl-MOM kinematics defined in the previous paragraph clearly have exceptional chan- 
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nels (e.g. pi + 



References [128 



p? — 0) and in this paper we generalize the non-exceptional Rl-SMOM schemes of 



30h to the four-quark operator. The choice of non-exceptional kinematics is not 



unique of course and in this paper we choose to study the Green function 

d{pi)s{~p2) ^ d{-pi)sip2) (12) 

with p\ = P2 = {pi — PiY = for a variety of momenta satisfying these conditions. In our 
notation below q = pi — pi- 

We briefly mention that we have previously investigated non-exceptional (or strictly speaking 
less exceptional) momenta for four-quark operators [280; here the operator was inserted only at 
a single point on the lattice and the method was less statistically precise than our current work. 
Chirality mixing in the four-quark operator basis arising in the infra-red p^ region was found to 



be strongly suppressed [|28|] . thus revealing the true, good chiral properties of DWF. However, the 
corresponding perturbative calculation to match this kinematic point to the continuum MS scheme 
was not available, and this was of largely academic interest in displaying the quality of Domain 
Wall Fermions. 

The remainder of the section is organised as follows. In the next subsection we introduce 4 RI- 
SMOM renormalization schemes, all of them defined with the kinematics of Equation (fT2l). In 
Subsection IIIIBI we calculate the perturbative matching factors relating ^vv+aa in the 4 RI- 
SMOM schemes with that in the MS-NDR renormalization scheme. We review some aspects 
of the non-perturbative renormalization of the lattice operator into a RI-SMOM renormalization 
scheme in Subsection IIII CI and finally in Subsection IIII Dl we combine the NPR computation and 
matching calculation to obtain the total renormalization factor relating the lattice and MS-NDR 
operators. 



A. RI-SMOM Renormalization Schemes for 



^VV+AA 



We follow the procedure which was defined for the renormalization of the four-quark operators in 



the RI-MOM Scheme 12711 . but now with the kinematics defined in Equation (1121) . We begin with 
the evaluation of the amputated four-point Green function A'^^'^g of the operator ^vv+aa, where 
05, /3, 7, and 5 are the spinor labels corresponding to the incoming s and d quarks and outgoing 
s and J quarks respectively and i, j, k, I are the corresponding colour labels. Analogously to the 
definition of the RI-MOM scheme, we impose conditions on the amputated Green functions at the 
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renormalization scale in such a way that they are automatically satisfied by the tree-level Green 
functions. To this end we introduce two projection operators P'^l^^^f^p yg, with X G {1,2}: 

^w!-p,rs = 256N{N+1) ^^fhairv)5r+if^)pa{rv7')sr] ^lAi . (13) 

P^UjS = 6V^(A^+1) [^^Pai^Sr+mpamSy] ^iAl , (14) 



where A'' = 3 is the number of colours. These projectors are constructed to give 1 when contracted 
with the tree-level result for ^'^'^yg given in Equation (l24l) below. 

In order to specify the renormalization condition on the operator we have to include a factor of 
^yz^ for every external quark line, where Zq is the wave function renormalization factor, and 
here again we use two possible definitions, called RI-SMOM and RI-SMOMy^, in Reference [|30|] . 
Here, we do not reproduce the explicit definitions in terms of the renormalization of the quark 
propagator, but note that they are chosen to satisfy the Ward Takahashi identities when combined 
with the renormalization conditions on the vertex function for the (conserved) vector current using 
two different projectors. Specifically in the SMOM-scheme 

^RI-SMOM ^_|_Tr[A(:^], (15) 

where the trace is over both colour and spinor indices, q is the momentum transfer at the vector 
current and Ay is the amputated two point function with the incoming (outgoing) quark having 
momentum pi (pz) with q = pi — pi and with p\= P2 = q^ chosen to be the renormalization 
scale. For the second scheme we use the same projector as in the definition of the RI-MOM 
scheme, but with the non-exceptional kinematics as above, 

Z,"*=^Tr[Ai;/']. (16) 

We label the renormalized four-quark operator in each of the four schemes by two labels {X,Y) 
with X = ^ or depending on which of the projectors Equation (fT3l) or (fT4l) are used for the 
vertex and similarly 7 = ^ or 7^ depending on which of the definitions Equations (fTSi) or (fT6l) are 
used for the wavefunction renormalization. Thus for example. 



where 
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We have introduced the subscripts R and B in Equations (fTTI) and (fTSl) to denote renormalized 
and Z^are (or lattice) quantities respectively. The remaining renormalized operators are defined 
similarly: 

4'-"' - (zr— , ' (19) 

^(l),ajS,75^fi,a/3,r5 

Z^^'^) = (zRI-SMOM)2 1 (20) 

G \ q ' pij,ki .ij,ki ^ ' 

^(2),aj8.y5^^S,a/3,75 

'^(2),a/3,y5^fi,a/3,75 
and in each case ^^^v;^^^ = Z^'^^ i^b, vv+aa, with X , 7 = ^ or 7^ . 

In addition to the four renormalization schemes defined above, we also use the standard RI-MOM 
scheme as the intermediate scheme in our conversion to MS. The reason for introducing sev- 
eral renormalization schemes is that it allows us some control over the lattice and perturbative 
uncertainties. After performing the perturbative matching to the NDR scheme, each of these in- 
termediate schemes should lead to the same value of the matrix element of i^yv^^- The spread 
of results obtained using the 5 schemes is therefore a measure of the uncertainties. In particular, 
since the matching coefficients from the intermediate schemes to the NDR scheme are currently 
available only at one-loop order (see Subsection IIIIB|) . the spread of results is an indication of 
the size of the higher-order terms. We now turn to the evaluation of the matching coefficient at 
one-loop order. 



B. Perturbative Conversion to the NDR Scheme 

In this subsection, we calculate the conversion (matching) factors between the four RI-SMOM 
schemes defined in Subsection IIII Al above and the naive dimensional reduction (NDR) scheme 
for the tsS — 1 operator = s^dsy^itd (where = y^{ \ —y') and we only consider the 

parity even component) using continuum perturbation theory at the one-loop level. The two-loop 
anomalous dimensions are also calculated to derive the renormalization group (RG) running of the 
operator in these schemes. 

We now use perturbation theory to convert the operators into the NDR schemes with the treatment 



of evanescent operators as in Reference [13 ih . as will be explained below. As explained above, for 



Bk the RI-SMOM schemes are defined in terms of projections of the amplitude d{p\)s{—p2) 
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P\ -pi -pi p\ p\ -pi -pi pi 











P2 ^ ^^Pl Pi ^ ^^-Pl Pi ^ ^^-Pl-Pl ^ ^^Pl 



FIG. 9: The four lowest order diagrams. Each circle represents the insertion of the current s'^d. The d or 
d {s or s) quarks have momenta itpi (±p2) and the flow of fermion number is denoted by the arrow. 



PiJ,P ^ ^-Pi,l,5 Pi,j,P ^ 





2 ^ 

{-P2,i,a^ ^\p2,k,Y p2,k,Y ^ ^\-p2,i,a 



FIG. 10: The lowest order diagrams, with spinor and colour labels exhibited. The notation is as in Figure|9] 

d{—p\)s{p2), where p\ = p\ = {pi — PiY = P^ with p\ ^ p2. For p^ in the perturbative regime 
there is no channel with soft momenta, thus reducing infrared effects. At tree level we have 
the 4 diagrams in Figure HI where the circles represent the two currents syj'^d, the arrows on the 
quark lines denote the flow of fermion number and the direction of the momenta are indicated 
explicitly below the corresponding momentum. Even though both momenta pi are ingoing and 
both momenta p2 are outgoing, it is convenient to introduce the minus signs and to think of the 
process as d{pi)s{—p2) d{—pi)s{p2) because then the signs also implicitly keep track of the 
spinor and colour labels (see Figure [TOl). Since the two currents commute, the first two diagrams 
are clearly equal as are the second two; thus we can rewrite the four diagrams in Figure|9]in terms 
of the two diagrams in Figure[lOl where the spinor (Greek letters) and colour (Latin letters) indices 
have now been indicated explicitly. The mathematical expression corresponding to the diagrams 
in Figure [To! is: 

2{(7l )a)8 (r/xL)y5 " (T^)ri8 ir^iL)ad (22) 

where the minus sign between the terms arises from fermion statistics. The Fierz identity (for the 
parity even component) 

(7L)a/3(rML)y5 = -(7f)r/3(7/iL)a5 (23) 
allows us to write the lowest order result as 



2(T^)„/3 {r^L)rS + 5a5kj} . (24) 
Writing the result in this way, the spinor structure is just that of the first of the four diagrams in 
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Figure m but the colour factor is different. It will be convenient in defining the projectors to take 
a trace in colour space, i.e. to multiply the expression in Equation (|24|) by 5,^5^/ and sum over the 
repeated indices. This gives a colour factor at lowest order of N^+N, where A'^ = 3 is the number 
of colours. 

We presented the above arguments explicitly because they generalize to the one-loop calculations 
below. Consider for example the 4 diagrams obtained by adding a gluon between the quarks with 
momenta labeled pi and —p2 in Figure [9l Each of these four diagrams can be Fierz-transformed 
into each other. It is therefore sufficient to calculate any one of the diagrams, but care needs to be 
taken in order to evaluate the colour factor correctly. 

Fierz identities are four dimensional relations whereas in NDR one works in Z) = 4 + 2e dimen- 
sions. This is the origin of the so called evanescent operators such as 



which vanish in 4-dimensions by the Fierz identity. Equation (1231) . Note the relative minus sign 
compared to Equation (l23l) due to the interchange of fermion fields. It is conventional to de- 
fine the NDR operators having subtracted the evanescent operators, i.e. using the 4-dimensional 
Fierz identities (analogously to subtracting the Euler constant and log(4;r) when defining the MS 
scheme). This is possible because the evanescent operators vanish in 4 dimensions and are there- 
fore proportional to £ and are only combined with the 1 / £ divergence. Their contribution is 



therefore independent of momenta. The evanescent operators are therefore removed by one 



counterterms, and must be included when evaluating the two-loop anomalous dimension yi 



In order to compare our result for the one-loop counterterms with Referenc e 113 ih we evaluate their 



oop 



3211 



coefficients. We use the same basis of three operators as in Reference plQ : in addition to Ei 
defined in Equation (l25l ) we introduce 

El = {s'r^rvrpPLd'){sJfffPLd')-{'^6 + 4£){s-'yj^d'){sJr^LdJ') (26) 
E3 = {s-^r^irvrpPLd^){s-YffPLd') -{16 + Ae){s-^y^d'){s-jY^,LdJ^ (27) 



where Pl= l-f ^M- In comparing our results with Reference Il3 111 the reader should note that 
we use D = 4 + 2e to denote the number of dimensions whereas the authors of Reference Isijl use 
D = 4-2e. 
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Pi 



-P2 





-Pl 



P2 



Pl 



-P2 




-Pl 



P2 



(al) 



(bl) 



FIG. 1 1 : The two independent one-loop Feynman diagrams to be evaluated. 



I. Evaluating the Diagrams 



There are two independent Feynman diagrams which have to be evaluated (see Figure [TTI) and we 
now present the results for these diagrams. The results are presented before taking the traces corre- 
sponding to the projection operators which define the RI-SMOM schemes, and so contain flavour 
and colour indices. The expressions for the remaining diagrams can then be readily obtained 
from those in Figure [H] by symmetries, except for the contribution of the evanescent operators to 
the one-loop counterterm which we also discuss later. Leaving the indices free also provides us 
with the flexibility to use a variety of renormalization schemes (such as the schemes defined in 
Subsection llll Al) which we exploit at the end of this Section. 
Diagram (al) gives the following result: 



16;r2 



16;r2 



5, 7 4/ 



l+2Co/^iT^ A 



+ 



2 hfR Pi+ hfR h 



®ypL 



3 



(28) 



log 



+ 



Co -4 



"1 ^ 2 ^ IP^ 



3 3 



(29) 



where Q = fi"{f) - Tr)^ ~ 2.34391 and ^(x) is the digamma-function ^(x) = T'{x) /T{x). In 
Equation (|28|) . X ®Y denotes X^^^Yy^, = /""(l + 7^) and t, is the gauge parameter defined so 
that = corresponds to the Landau gauge and = 1 to the Feynman gauge. It will prove to be 
a convenient shorthand to define A^^jg yg as in Equation (|29l ). 
The expression for diagram (bl) is 



(l-^)Tf®rpL 



1 ^_2(l-log2) 
4^Ai' 3 

, 4(l-log2) 



+ 



+ 



(30) 
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-P2 



Pi 



P2 





(al) 





(a3) 



-Pi 



P2 



-P\ 



-P2 



Pi 



-P2 



P\ 



P2 




(a2) 




(a4) 



-Pi 



P2 



-Pi 



-P2 



FIG. 12: Four one-loop diagrams whose Feynman integrals are given by that of diagram (al) in Figure [TTl 



Af^YpLAYn / l + 81og2 _ 



4iog2- r 



Diagrams (al) and (bl) in Figure [H] are not the only ones which need to be evaluated but, apart 
from the subtlety associated with the evanescent operators (which we neglect for the moment but 
to which we return shortly), they are the only ones for which the Feynman integrals need to be 
evaluated. Consider first the four diagrams in Figure [l2l in which one end of the gluon is attached 
to the quark labeled with momentum p\ and the other to one with momentum ±P2- The results of 
the four diagrams in Figure [12] can then be deduced by inspection: 

{al)=Aa^^ysCF8ijdki\ (al) = -Ayp^asT-jTg; ^^^^ 

(a3) = -Ay^ asCF^iiSkf, (fl4) =Aap y§T"jTg . 

To these must be added the contributions from the four diagrams in which one end of the gluon is 
attached to the quark with momentum —pi. These are obtained from the results in Equation (|3T1) 
by making the substitutions a ^ 7,/3 -v^ 5,i -v^ k, j -v^ I, and the sum of the eight diagrams is to 
be multiplied by 2 to include the diagrams obtained by interchanging the two currents. In this way 
we obtain a total answer for the 16 diagrams in which a gluon is attached to quarks of different 
flavour 



Ca = 2{Aapj8+Ay5,ap){CFSijdki + T^"T^j) - 



(32) 



1 



\6n^ e 



7 tree 



1 . 
—1 

N 



7 tree 



tree 



The last term contains the contribution from the evanescent operators which we have ignored up 
to now in this discussion. They arise because in rewriting the divergent terms in terms of the 
spinor structure (^f )ajs(7pL)y5 or {'yL)ap{YpL)'YS we have used the spinor Fierz identities which 
are not valid in D = 4 + 2£ dimensions. These contributions only arise in the presence of the e 
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Pi ^ -Pi Pi , -Pi 





-P2 ^\ P2 Pi -Pl 

^ (bl) ^ ^ (b2) ^ 

Pl -Pi Pl -Pi 





-P2 ,,->-SVBTriririrB-!nfS^ P2 P2 ,,-->-S=BWinrrin5-!f8-^ -P2 

(b3) (b4) 

FIG. 13: Four one-loop diagrams whose Feynman integrals are related to that of diagram (bl) in Figure [TT] 

ultraviolet divergence and are hence straightforward to identify. When evaluating the conversion 
factor between the RI-SMOM and NDR schemes, we will use projection operators which have 
some symmetry in the indices and which effectively simplify the expression in Equation (|32l) . 
Next we consider the 8 diagrams whose Feynman integral is given by the expression in Equa- 
tion ( |30l ). Four of these are shown in Figure [13] and the remaining 4 are obtained by switching the 
two currents (and are equal to those in Figure [T3l). The result for each of the diagrams (b2)-(b4) 
can be deduced by inspection from that for (bl) given in Equation (|30l) and for the total contribu- 
tion from the 8 diagrams we find: 

g2 A^_i /i ,n4-aw^' 4(l-log2) 

^KprSujki f 1 + 8 log 2 , ^ 4 log 2 - 1 



where 



^it^^^i^^-^^-^)^^^ (33) 



2. Tlw Conversion Factor 

Having kept the external colour and spinor indices uncontracted in Subsection llll B 11 we are in a 
position to determine the conversion factors relating the A5 = 2 four-quark operator defined in the 
four RI-SMOM schemes to that in the NRD scheme. The conversion factors, Co % are defined 
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by 



^vv^AA(M) = C^b/\p'/^^') 4v+AAiPl (35) 
where for convenience at this stage we keep p as the renormalization scale in the RI-SMOM(X, Y) 
schemes and jl as the renormalization scale in the NDR scheme. Since in this subsection we 
are only concerned with renormalized quantities we drop the subscript R denoting renormalized. 
From the definition of the RI-SMOM renormalization schemes given in Equations (fTSi) - (|2T]) we 
see that the conversion factors can be obtained from the equations 



^ / pi.i,kl 



A 



NDR, ij,kl 



c 



Bk 



(36) 



where, as throughout this paper, A represents the amputated Green function. Q are the conver- 
sion factors relating the wave-function renormalization factors in the MS scheme and that in the 
RI-SMOM scheme labeled by Y, dp = Zf^/Z^. At one-loop order these were already obtained 
in Reference yOfl, 



RI-SMOM 



C 



RI-SMOM^^ 



1 + 



1 + 



16;r2 



1 



1 



+ 



Cf 



1 



21og^ 



-Co 



+ 



(37) 
(38) 



where Cf denotes the Casimir operator in the fundamental representation of SU(A'^). These results 



have recently been extended to two loops 1133 . 



We now sketch the calculation of the conversion factor for the RI-SM0M(7^, ^) scheme and then 
present the results for the other three RI-SMOM schemes. The renormalization condition in Equa- 
tion (|36l) with the projector of Equation (fT3l) in the (y^i,^ scheme can therefore be written in the 
form 



^^RI-SMOM N 2 pij,kl . NDR, ij,kl 



c 



Bk ■ 



(39) 



non-except. 

From Equation (|39l) . together with the expressions in Equation (|32l) , (|33l) and (|37] ) we can evaluate 
the conversion factor between the (7^,^ and the NDR scheme. 
There are 3 contributions to the conversion factor: 

1. The total contribution from diagrams such as those in Figure [T2l above, in which the gluon 
is exchanged between a strange quark or antiquark and a down quark or antiquark, is: 



g2 (A^_i)(A^ + 2) 
16;r2 



log - 1 + I Oyyl^ip) 
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1 



+ 



-(S + 2^)Ei-^E2 + ^E3 



I6n^ 2e 

where A'^ = 3 is the number of colours ((3 — Cq)/2 ~ 0.328046). 



(40) 



2. The corresponding contribution from diagrams, such as those in Figure [13] above, in which 
a gluon is exchanged between quarks of the same flavour (i.e. the two strange quarks or the 
two down quarks), is: 



' ^ ^ '(3 + ^)log^ + 121og2-7 + 2^(21og2-l)|oiJ;fi,(;.) 



I6n^ N 

^2/1 1 



16;r2 V 4e 



1 



(1-^)^1 • 



(41) 



3. Finally we have the contribution from the quark wave-function renormalization: 



(42) 



Before presenting the final result we make two observations: 



1 . The total term with evanescent operators is 

-2 1 /I 1 ^ 

.{E2,--E2)-5Ex 



\6n^ £ V2 



(43) 



This term is eliminated by introducing counterterms which are equal and opposite to this. 
The result agrees with (2.15) and (2.22) of Reference [31] (recall again that we are using 

n 

D = 4 + 2£ and the authors of [31] are using D = 4 - 2£). 



2. The total logarithmic term is 



(44) 



which agrees with the known anomalous dimension. 
The final result for the conversion factor C^^'*^^ is given by 

„2 



c 



{7)14) 



Bk 



1 + 



16;r2 



\- (9-31og^-121og2 ) -8 + 121og2 + 31og^-A^ 



+ ^(i(Co-41og2)-l-^ + 41og2 + ^(l-Co) 



+ 



N=3 



1 + 



21og^ + 81og2-8 + ^ (^l_£co + ^log2 



+ 
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1 + 



16;r2 



21og^- 2.45482- a 1.05812 



+ 



(45) 



The remaining three conversion factors are obtained from equations Equations (|32l) . (|33] ) and dTT]) 
or (|38] ) in a similar way and we only present the final results. For the (7^, 7^) scheme we find 



C 



Bk 



1 + 



16;r2 



1 



8-121og2-31og^ 
11^ 



8 + 121og2 + 31ogi-y 



^(1 



,2^..+C„-81og2)-i-^ + 41og2 



+ 



iV=3 



1 + 



16;r2 



16 



21og^ + 81og2- — + Co--log2 



1 1 



8 



3 3 



+ 



1 + 



16;r2 



21og ^ + 0.21 1844 + 1, 0.733757 



+ 



For the remaining two schemes we use the second projector in Equation (fT4l) and impose 



(C 



pijM ^NDR, ijM 



C 



Bk 



non-except. 



again with q = pi— pi and Pi = P2 = = p^ ■ The conversion factors are 



and 



1 + 



16;r2 



1 

N 



9-31og^-121og2 +121og2-9 + 31og 



+ ( -(Co-41og2)-Co + 41og2 



+ 



N=3 



1 + 



1 + 



16;r2 



21og^ + 81og2-6 + <^ ( -Iog2--Co 



8 



21og ^ - 0.454823 + ^ 0.285788 



+ 



+ 



C 



Bk 



1 + 



16;r2 



1 



8-121og2-31ogi^ +121og2-9 + 31og^+A^ 



+ ^ ( ^(l+Co-81og2)-Co + 41og2 + ^(Co-i; 



N=3 



1 + 



1 + 



16;r2 



10 



21og^ + 81og2-y + <^ ( -Iog2 + -Co-- 



8 



21og ^ + 2.21 1844 + ^ 2.077664 



+ 



(46) 



(47) 



(48) 



(49) 



The results for the four conversion factors for the RI-SMOM schemes together with that for RI- 
MOM are summarized in Table lVIII 
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Scheme for 
four quark operator 


Cb, for ^ = 


RI-MOM 


1 + ^(0.87851...) + ^(a,^) 


(7m 


1 + ^(-2.45482...) + ^(a,^) 




1 + ^(0.21 184...) + ^(a2) 




1 + ^(-0.45482...) + ^(a2) 




1 + ^(2.21 184...) + ^(a2) 



TABLE VII: Summary of the conversion factors (in the Landau gauge) of the four quark operator from the 
RI-(S)MOM schemes to the MS[NDR] scheme. 

3. Two-Loop Anomalous Dimension 



We follow the conventions of Reference fi3 lil and define the anomalous dimension y of the renor- 
malized operator O by 

dn(n\ 

(50) 



M^--r(M)W, 



where /i is the renormalization scale. Expanding 7 as a perturbation series 



r(Ai) 



(16;r2) 



\6k^ ) 



(51) 



the one and two-loop coefficients in the MS-NDR scheme (called NDR in the following) are 113 511 

y(0)NDR ^ 6_ 6 A'=3 4 



r 



1)NDR 



22 57 39 19,, /2 2 \ Af=3 ^ 4 

N + Hf = -7 + -nf, 

3 2N^ N 6 3nJ 9 ^' 



(52) 
(53) 



where ny^ = 3 is the number of flavours contributing to the running in the region of interest. 
Now let the conversion factor between the NDR scheme and a scheme A which is defined in the 
Landau gauge so that the gauge parameter is not renormalized be given by 



l + ^^ArA^NDR + - 



(54) 



In the following we consider for the 5 schemes A G {RI-MOM, (7^1,^), (7^1, 7;x), {^,^), (^, 7/i)}- 
From Equation (l45l) we see that ArRi.sMOM^NDR — —2.45482 and from Section 5 of Refer- 
ence [,3 1.1 we read 



ArRi.MOM^NDR = -7 + + 12 M - ^ j log2 V 0.87851 1 . 



(55) 
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For the one-loop anomalous dimensions the equation y^^^^ = y(0)'^DR j^glds and the relations 
between the two-loop anomalous dimensions are given by 



Jl)A ,,(1)NDR o« A»- 
7=7 - ^PoArA^NDR , 



(56) 



where j8o is the one-loop coefficient of the QCD j8 -function which is defined by 



i8 



-/3o 



An 



/3i 



An 



(57) 



with 



/3o 
/3i 



11 



-N--n 



34, ,2 13. 



N\nf.^ 



(58) 
(59) 



and 0(j(/i) = g^{iJ.)/{An) is the strong coupling constant. In this way we obtain in the Landau 
gauge 



7' 



1)NDR 



Y 



,(l)RI-MOM 



57 39 22 19 2 
57 39 176 



1 

1 

N 
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r 



.(1)(7m, 



2N^ N 
1 /26 

57 39 220 



^ +881og2-FA^l ^ 



161og2 ) - y + 161og2 



A'=3 _17 

"/=3 3 
-88 log 2 

A'=3 



/!/ = 3 



-21.4799, 



r 



(i)(rM.rp) 



1 /34 

57 39 
+ ■ 



111 \ 22 o 

^ +881og2 + A^( — -881og2j+— 



yv=3 

~ 38.5201, 

«/=3 



161og2 ) -10+161og2--A^ 



2A^2 ^ 



66 + 881og2 + iV 



111 



-88 log 2 



^ (10 - 161og2) + 161og2 - 10 



^=3 

~ -9.47986, 

n/'=3 



7' 



1/34 \ 34 



A'=3 

~ 2.52014, 

"/=3 



■(i)(^,rM) 



57 39 /377 \ 22 . 

^ + --66 + 881og2 + ^(— -881og2 --A^2 



i (10 - 161og2) - ^ + 161og2 + 



N=3 

- -45.4799. 

n/=3 



(60) 



(61) 



(62) 



(63) 



(64) 



(65) 
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In Reference [|32l 13611 a factor has been introduced to convert the results to the renormalization 



group independent (scale invariant) value defined by 



4f(«/) = £0i^(Ai,«/)zi(Ai,'^/): 



(66) 



where A again labels the scheme. At next-to-leading order the contribution to the evolution of the 



operator is written in terms of a quantity called 7^ 



("/) 



-r'"V(2/3o) 



as defined in Appendix D of Reference fesll . In the notation used here it is given by 



With = 3 we find in the Landau gauge 



2Po 2/32 



r 



.(3) 
'NDR 

(3) 

RI-MOM 

r(3) 

(3) 



13095 -1626n/ + 8n^ 

6(2n/-33)2 njts 
17397- 2070/1/ +104n2 



1.89506. 



6(2n/-33)2 
39177 -471 On/ +184^2^ 



/ 



J 



6(2«/-33)2 
7251-866n/ + 40n^ 



+ 81og2 ~ 2.77357, 

nf=3 



+ 81og2 ~ -0.55976, 

nf=3 



2(2„,-33)^ +81og2_^c.^ 2.10691 



(3) 



26109-3126/1/ + 136^2, 



7 



7 



(3) 



6(2/1/ -33)2 
2895-338/z/ + 24/i2 



+ 81og2 ~ 1.44024, 

«/=3 



0^0 oo^2 +81og2 ^ 4.10691. 
2(2/z/-33)2 ° „y=3 



(67) 



(68) 



(69) 
(70) 
(71) 
(72) 
(73) 
(74) 



The first two results in Equations ( l69l) and (iTOl) can be taken from Reference [l32ll and agree with 
(D4) and (D3) respectively in Reference [i28ll . 



C. Volume averaged vertex functions 



In contrast to earlier RBC-UKQCD publications [1280, in the present study we have developed 
volume-source NPR for four quark operators with a generalised momentum configuration. As 
will be demonstrated below, this volume averaging greatly improves the statistical precision. The 
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technique is similar in style to previous analyses introduced for bilinear operators by the QCDSF 
collaboration [|37D. The advantage of the method arises from the fact that the amputated vertex 
functions are evaluated with the operator insertion averaged over all lattice sites, as opposed to 
the single-point source operator insertion. The resulting statistical errors are tiny and systematic 
effects like Ga, breaking lattice artefacts dominate. These must be included in the error analysis or 



removed using, for example, the techniques of [13 811 (which we also do in this study). 



We define the four momentum source, used on a Landau gauge-fixed configuration, as 

7]pW=e'>''-'5,,5«^, (75) 



where z, j and a, /3 are color and spinor labels respectively and the momenta take the values 

1% 



(76) 



where n is a four-vector of integers. 

On a given gauge field {x) we solve the equation 

M{x,y)Gp{y)=np{x\ {11) 

and M is the domain wall fermion matrix with (5 — M5) 1 on the site diagonal portion. 
In performing the NPR, as explained above, we select two momenta p\ and p2 satisfying p\ — 
P2 = {p\^ PiY- In order to reduce the artefacts arising from the breaking of ^4 symmetry, we 
selected values for Pi= P2 = {pi — Pi)^, such that while still satisfying the Fourier constraints 
we best minimise Y.iP^ as documented in Table IVIIIl Alternatively, following ref, 113 811 . we may 



impose twisted boundary conditions [|39l - l44ll on the quark fields 



^(x + L) =e'%(x) where B^ = ^ (78) 

Equation (1771) is then modified to 

M{x,y)Gp{y) = r)p{x) where G{y,p) ^ e-'^'yCp+Biy) (79) 

Thus by varying the twist angle 6 we can vary the magnitude of the momentum without changing 
the direction. Our choices of p and B are documented in Table |IXl The particular choices here 
are the non-exceptional directions that minimise Y^iPf- We choose the components of B equal and 
always in the same direction as p: for example if p = (0, 1, 1,0) then B = £(0, 0, 0,0) . 
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24^ X 64 pi 


Pi 


32^ X 64 


Pi 


(0,4,4,0) 


(4,0,4,0) 


(3,2,2,2) 


(3,2,-1,-4) 


(1,2,2,8) (-2,-1,2,8) 


(4,2,2,0) 


(4,0,-2,4) 


(1,4,2,8) 


(2,-1,4,8) 


(4,4,3,2) 


(4,3,-1,-8) 


(2,2,4,0) 


(4,-2,2,0) 


(4,-5,0,-6) 


(4,0,-5,-6) 


(2,3,2,8) 


(3,-2,2,8) 


(-4,-1,-4,2) (-4,-4,1,2) 


(-3,1,1,8) 


(1,1,3,8) 







TABLE VIII: Non-exceptional discrete momenta used for the evaluation of amputated Green's functions 
in our NPR analysis. The momenta here are listed in {x,y,z,t) order for our 24^ x 64 and 32^ x 64 lattices. 
The integer Fourier mode numbers {n,} are given and the lattice momenta are related via api = The 
exceptional momenta used correspond to p2 = p\ for the same set of momenta. 



24^ X 64 


P\ Pi 


e 




(-3,0,3,0) (0,3,3,0) 
(-4,0,4,0) (0,4,4,0) 


^n:n = {-2,1..., 12} 

3 
2 


32^ X 64 


P\ Pi 







(-3,0,3,0) (0,3,3,0) 
(-4,0,4,0) (0,4,4,0) 
(-5,0,5,0) (0,5,5,0) 


1 

4 

3 3 

4 ' 8 

5 3 
8 ' 8 



TABLE IX: Non-exceptional momenta and twist angles used for the evaluation of amputated twisted 
Green's functions in our NPR analysis. The momenta here are listed in (;c,y,z,f) order for our 24^ x 64 and 
32^ X 64 lattices. The integer Fourier mode numbers {«,} are related to the lattice momenta via api = 
The momentum added by the twist, B, is determined by the twist angle B giving apt = j^e 
exceptional momenta used correspond io p2 = p\ for the same set of momenta. 

We now form phased propagators 

G'p{x) = We-'>^ = Y^M-\x,y)e'P<y~''^ . (80) 

y 

With twisted boundary conditions this equation is generalized to 

y 
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so that the phases are properly accounted for and the following discussion holds for both twisted 
or untwisted propagators. For each configuration we form unamputated bilinear and four quark 
vertex functions for generic Dirac structure F: 



(82) 

-I ij,ap 



and 



E(r5(G;,w)tr5FG;,w)^.^.^^(r5(G;,W)V5rG;,w)^^^^^. (83) 

Here, external colour and spin indices are left free for later amputation. We use the kinematics ex- 
plained in Section lTlIBl in which the four-point functions have two legs with incoming momentum 
pi and two with outgoing momentum p2. 

A single 12x12 object is written out for each configuration and momentum point for the bilinear 
vertex functions, and al2xl2xl2xl2 object for the four quark operator. For convenience, we 
use a single 12 valued index below to represent both color and spin. These building blocks enable 
the accumulation of the following ensemble averages 



5).. = UiG'Mab)^ (84) 

{Vr{puP2))ab = (£{Y5{G'p,)Hx)r5rG'p^{x))j, (85) 

wripuPi) = (i:(75(G;jnx)75rG;,^w)^^^(r5(G;jnx)r5rG;^w)^;. (86) 

These ensemble averages are then used to construct the amputated vertex functions for bilinears 

^Minear ^ 75 (G^ ) " Vs VfIt^I , P2) (G;^- ^ , (87) 

where F e {A, V, 5, P, T} and for four quark operators 

4' = (r5(G;j-V5)^^, (r5(G;,)-V5)^^wr'(pi,;'2)(G;,),-,.i(G;j,-:,i (88) 

v/here r e {VV ±AA,SS±PP,TT} . 

Finally the A^'^ are contracted with the projectors defined in Equations (fT3l) and (fT4l) . 



D. Lattice Results for the Renormalization of Bfc 



While the methods summarized in the previous section can be directly applied to the case at hand, it 
is important to adopt a strategy which depends on amplitudes which can be accurately determined. 
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For example, it is useful to directly calculate the ratio of renormalization factors in the scheme 
S, Z^^^^^/Z^, which is needed for the ratio of the four quark matrix element to which enters 
the actual definition of Bk because the common factor of Z^ appearing in the lattice calculation of 
^^vv+AA cancels in this ratio. (Here Zq is renormalization factor for the domain wall quark 

field which is central to the RI-MOM approach but may introduce large systematic errors if it is 
identified as the coefficient of a momentum-dependent term in the lattice quark propagator.) 
Thus, we transform our lattice-normalized result for to one normalized in the scheme S by 
multiplying by the ratio 

Zl = ^-^=(jtj]' , (89) 



\ CTyv+AA / m^O 

where T^^y^^ is the projection of the amputated Green function, A^'^, with a projector from Equa- 
tions (fT3]) and (fT4l) corresponding to the renormalization scheme S, and Ty = ^ is the appropriate 
projection of the amputated vertex function of the local vector current Ay. Here either the local 
vector or axial current can be used since their difference is expected to be of order m^g^. 
We compute Zg^ in each scheme using Equation (l89l) . The twisted momenta are given in Table HXl 
For the 1 ensembles the lattice momenta approximately span the physical range 4.0 GeV < p < 
ll.OGeV^. On the 2 ensembles the momenta span 3.25 GeV^ < p~ < 9.0 GeV^. The overlap 
region, 4.0 GeV^ < < 9.0 GeV^, will be used for continuum extrapolations. 
We perform a linear extrapolation of the results to the massless limit using data with quark masses 
corresponding to the dynamical light-quark masses m;. We do not observe any statistically rele- 
vant mass dependence in Zg^^. Since we are restricted to a single sea strange quark mass in our 
computation, we cannot perform a chiral extrapolation for the third active flavour. This mismatch 
between the mass -independent renormalization schemes and the finite sea strange quark mass is 
included in our error budget. 

The lattice data in the chiral limit is converted to the NDR scheme at the renormalization scale 
jl = 2 GeV or /i = 3 GeV using the perturbative results from Section III B. 

Several additional inputs are required: we define the three flavor coupling a, from the PDG 2010 
central values a,(Mz) = 0.1184(7), mp = 4.19+J^ GeV and mp = 1.21^1 GeV by using the 
four-loop running down to our renormalization scale and matching across flavor thresholds. We 
combine this four-loop and 2-1-1 flavour with the two-loop anomalous dimensions to obtain the 
Wilson coefficients for both scheme change to MS, and to obtain the 2-1-1 flavour RGI operator. 
The perturbative contribution to the momentum scale dependence is divided out, and the data for 
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is displayed in Figure [14] and [151 The remaining dependence is a source of systematic error 



and is discussed in detail in Section UlID 1[ 
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FIG. 14: We can use the the perturbative running to convert the chiral limit of the ratio (1891 ) to MS at 2 GeV 
for each using Zg^{p^) x "'^^g(^2~f^''y,'^."^3y^^ ■ This is displayed for all five intermediate MOM schemes 
5 on the 2 ensemble set (24^, a^^ = 1.73, GeV lattice). The top two panels correspond to the original 
RI-MOM as the intermediate scheme and the other four rows coiTcspond to the schemes of Section HUB I 
The left-hand panels show the data with the momenta of Table lVIlTl and the right-hand panels show the data 
using the momenta in Table IIXI accessible with the use of twisted boundary conditions. The scatter due to 
the 0(4) symmetry breaking in the left hand panels is absent in the right-hand panels. For this reason we 
use the data with twisted boundary conditions for our analysis. 



1. Systematic errors due to renonnalization 

In Tables |X| , [Xlj and IXIII , IXIIII we summarize the results and the error budget for the schemes 
described in Section ITIIAI There are six main contributions to the total error 

1. Statistical errors. These are denoted by the label "stat" in Tables IX] - tXIII[ 

2. Errors due to the breaking of ^4 symmetry. As explained below we eliminate these errors 
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FIG. 15: We can use the the perturbative running to convert the chiral hmit of the ratio (1891 ) to MS at 2 GeV 
for each using Zg^{p^) x ^'^^^(^£^^''^^-3^^^ ■ This is displayed for all five intermediate MOM schemes 
on the 1 ensemble set (32^, a^^ = 2.28 GeV lattice). The top two panels correspond to the original RI- 
MOM as the intermediate scheme and the other four rows correspond to the schemes of Section ITlIB I The 
left-hand panels show the data with the momenta of Table IVIIII and the right-hand panels show the data 
using the momenta in Table IIXI accessible with the use of twisted boundary conditions. The scatter due to 
the breaking of 0(4) symmetry is smaller on this finer lattice. 



by evaluating the Green functions using momenta which are made accessible by the imple- 
mentation of twisted boundary conditions. These are therefore absent in Tables |X] - IXIII[ 

3. Uncertainty in the values of the lattice spacing. We denote these by a^^ in Tables IXl - IXIIIi 

4. Uncertainties due to infrared chiral symmetry breaking effects. These are only significant in 
the RI-MOM scheme where one manifestation is the difference in the values of Ay and A^. 
We therefore label these effects by V - A in Tables ISSni 

5. Errors due to the fixed sea strange-quark mass when defining mass-independent renormal- 
ization schemes. We label this by nis in Tables |X] - IXIII[ 

6. Error due to the truncation of the perturbation series in the matching. We label this by PT. 
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scheme 


MOM 


SMOM (Yf^Yn) 


SMOM (7^,^) 


SMOM (^,7^) 


SMOM (^,^) 


Z^DR(2GeV) 


0.95541 


0.96089 


1.03838 


0.92164 


1.00028 


Stat 


0.00151 


0.00046 


0.00093 


0.00104 


0.00036 


a-' 


0.00045 


0.00052 


0.00211 


0.00030 


0.00129 


Ms 


0.00846 


0.00221 


0.00386 


0.00174 


0.00151 


V-A 


0.00551 


0.00014 


0.00013 


0.00010 


0.00014 


Total 


0.01022 


0.00232 


0.00450 


0.00205 


0.00202 



TABLE X: EiTor budget, without the perturbative truncation (PT) eiTor, for Zg^ (2 GeV) on the 1 ensem- 
ble set (j8 = 2.25, 32^ lattices.) 



scheme 


MOM 


SMOM (7^,7^) 


SMOM (7^,^) 


SMOM (^,7^) 


SMOM (^,^) 


ZNDR(3GeV) 


0.93453 


0.94284 


0.99252 


0.91681 


0.96698 


Stat 


0.00030 


0.00017 


0.00034 


0.00038 


0.00013 


a-i 


0.00058 


0.00049 


0.00137 


0.00004 


0.00086 


nis 


0.00181 


0.00048 


0.00039 


0.00024 


0.00009 


V-A 


0.00188 


0.00002 


0.00002 


0.00002 


0.00002 


Total 


0.00269 


0.00070 


0.00147 


0.00046 


0.00088 



TABLE XI: Error budget without PT error for Z^^^i3GeV) at j8 = 2.25 (32^ lattices) 



Since we estimate this error by comparing the results obtained in different schemes, it is 
absent in Tables |X] - IXIIII where errors in individual schemes are presented separately. 

We define the central value for Zg^ through a linear interpolation in (a/?)^ to the same physical 
scale = /i^ on both ensemble sets, and this is our chosen MS renormalization scale /i. We take 
the continuum limit of the renormalized matrix element, removing the lattice artefacts. This ap- 
proach differs from earlier work in our collaboration i280 where the values of the renormalization 
constants extrapolated to = were used. 
We now consider the sources of systematic error in more detail: 
1^4 breaking: 

The use of volume sources leads to tiny statistical errors and as a result the scatter of the points 
around a smooth curve in {ap)^ becomes a prominent source of uncertainty. This is illustrated by 
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scheme 


MOM 


SMOM (7^,7p) 


SMOM (7^,^) 


SMOM (^,7^) 


SMOM (^,^) 


Z^DR(2GeV) 


0.92578 


0.93731 


1.01350 


0.89936 


0.97621 


Stat 


0.00028 


0.00010 


0.00032 


0.00027 


0.00011 


a-' 


0.00049 


0.00064 


0.00225 


0.00013 


0.00140 


Ms 


0.00757 


0.00393 


0.00445 


0.00054 


0.00180 


V-A 


0.00750 


0.00021 


0.00026 


0.00021 


0.00026 


Total 


0.01067 


0.00399 


0.00500 


0.00065 


0.00230 



TABLE XII: Error budget without PT eiTor for Z^^^ {2GeV) on the 2 ensemble set (j8 = 2. 13, 24^ lattices). 



scheme 


MOM 


SMOM (7^,7^) 


SM0M(7p,^) 


SMOM (^,7^) 


SMOM (^,^) 


ZNDR(GeV) 


0.90444 


0.91983 


0.97455 


0.89147 


0.94672 


Stat 


0.00066 


0.00010 


0.00029 


0.00027 


0.00011 




0.00076 


0.00051 


0.00131 


0.00007 


0.00084 


Ms 


0.00347 


0.00181 


0.00164 


0.00148 


0.00063 


V-A 


0.00203 


0.00003 


0.00012 


0.00009 


0.00012 


Total 


0.00415 


0.00188 


0.00213 


0.00151 


0.00106 



TABLE Xm: Error budget without PT error for Z^^^{3GeV) on the 2 ensemble set (j8 = 2.13, 24 
lattices). 



a comparison of the left and right-hand plots of Figures [14] and \T5\ The scatter in the left-hand 
plots, which correspond to Fourier momenta given in Table IVIIIl can be attributed to artefacts 
which appear due to the breaking of rotational symmetries on the lattice. In previous studies they 
have been hidden due to the statistical noise and the averaging over all degenerate p^. In a recent 
paper 113 811 it has been shown how this scatter can be avoided using twisted boundary conditions. 
Instead of using the Fourier modes, we introduce twisted boundary conditions and use momenta 
which are equivalent under the hypercubic group on each lattice spacing. This eliminates the 
spread due to the breaking of ^4 invariance. This expectation is confirmed in the right-hand plots 
in Figures [14] and [T5l where we use the twisting angles specified in Table ]IX] and we therefore use 
the twisted data exclusively in this analysis. Of course, the 0{a^) errors still remain - we have 
simply chosen a single orientation for the lattice momentum. The twisting allows us to deal with 
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these discretisation errors by taking the continuum limit of a fixed observable with a controlled 
Symanzik expansion. 
Uncertainty in the lattice spacing: 

In order to obtain the renormalization constants at a given physical scale we use our measured 
values of the lattice spacings = 1.73(3) and agj^ = 2.28(3) [19j. The central values quoted 
above for the renormalization constants are obtained using the central values for and the errors 
are estimated by recalculating Zg^ using a^^ + Aa^\ where Aa~^ is the error in the inverse lattice 
spacing, and taking the difference for the estimated uncertainty. 
Infrared chiral symmetry breaking effects: 

In the original RI-MOM scheme the difference between the bilinear vertex functions of the vector 
and the axial vector current is significant [28]. We perform separate analyses using or ^(Av + 
Aa) in Zg^, as these differ for the original RI-MOM kinematics due to infrared chiral symmetry 
effects. We include the difference as a systematic error and take the ratio with Ay as the central 
value. This was estimated to be one of the largest sources error in our previous Rl-mom work, 
but we now find that there is no measurable difference between the two cases for the new SMOM 
schemes, 
m,: 

We associate an error due to our treatment of data with sea strange quarks near their physical 
mass while using a mass-independent scheme when converting to MS. This can be estimated 
by measuring the slope of the data with respect to the simulated light-quark masses in the chiral 
extrapolation of vertex functions. We take one half of this slope, as there is now a single flavour, 
and multiply by the simulated strange quark mass to obtain the systematic error. This error is 
rather small for the non-exceptional momentum schemes which have a mild mass dependence. 
Perturbative truncation: 

For each scheme a perturbative truncation error arises because we only know the perturbative 
running to some fixed order. Estimating this error is necessarily subjective as a rigorous estimate 
would require us to know the unknown higher order terms. 

At fixed order there are two possible approaches that may be advocated as being reasonable es- 
timates of this error. Firstly, notional convergence of the perturbative series could allow one to 
estimate the error as either the last term in the series, or perhaps a", where n is the order of the 
first unknown term, or even (^)" according to subjective taste. These differ greatly, however for 
our preferred scheme SMOM(^, ^) the last term is around 0.8%. 
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Another approach is to compare the results obtained using different schemes to the order at which 
we know the results, and consider that any discrepancies between the schemes after the well- 
controlled continuum limit has been taken are indicative of the residual perturbative uncertainty. 
Here again some subjectivity enters through an assessment of which and how many schemes 
should be considered, however this is a promising approach which we adopt. 



perturbative running than the other schemes. Here we also find that the residual dependence 



our ensembles with a larger volume. This indicates that in the continuum limit, the SMOM(^,^) 
scheme is best described by the perturbative running, and we take the result in this scheme as 

(3) (3) 

our central value. We note that of our schemes J]^^^-^ was closest to ^^dR' ^'^^ ^^^^ therefore 
consistent with the small size of the perturbative correction needed to change scheme. For the 
error, we take the difference between the two schemes that are best described by perturbation 
theory in Section IfflDBi namely the difference between the SMOM(^,^) and S MOM (7^, 7^) 
schemes. 

We examined alternate strategies involving a weighted average of the results in all the schemes. 
This selects the schemes best described by perturbation theory, and deweights those poorly de- 
scribed by perturbation theory. Here the relative weight might be determined by the slope of each 
scheme after removing perturbative running. We find that in this case the overall error is slightly 
smaller than that obtained from the difference of the results in the SMOM(^, ^) and SM0M(7^, 7^) 
schemes, and so we adopt the latter as the more conservative error. 

We also note from our tables that at the higher scale the difference between schemes is smaller. 
For example on our finer lattice, i.e. closer to the continuum limit, we find that the rms error 
between the different schemes is reduced from around 0.04 to 0.03 as we go from 2 to 3 GeV. At 
a sufficiently high scale and in the continuum limit all schemes should give the same result. Since 
the difference between schemes is a major systematic error and we believe we have good control 
over lattice artefacts by taking the continuum limit, we prefer to compute Zg^ at the higher scale 
of 3 GeV. The non-perturbative conversion factor to go from 2 to 3 GeV in a variety of schemes 
will be presented in a later section. 

Finally, as a result of using a formulation of lattice QCD with good chiral properties we have no 
systematic error associated with operator mixing, as we explicitly demonstrate in the following 
subsection. 




scheme was better described by two-loop 



for the SMOM(^,^) scheme is the smallest, and in Section ITlI D 31 confirm the analysis of 113 811 on 
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2. Operator mixing 

The four-femiion operator Ovv^aa renormalizes multiplicatively when chiral symmetry is pre- 
served. This holds, for example, for lattice regularizations which preserve chiral symmetry and 
mass-independent renormalization schemes. In Reference [Zsl it was shown that the original RI- 
MOM procedure, with four identical momenta in the four-point vertex function, does not lead to 
vanishing mixing with the remaining elements of the basis of dimension six operators. Already in 



Reference 112811 it was pointed out that schemes with non-exceptional momentum configurations 
Pi= Pi = (Pi ^PiY give mixings consistent with zero. The application of momentum sources to 
this problem dramatically decreases the statistical error on the mixing coefficients. Therefore we 
are able to give more stringent bounds on the residual mixing which is expected to be of 0{am^g^) 
for domain wall fermions. In Figure [16] we present results for the mixing coefficient Zvv+aa.x, 
where X = VV -AA,SS-PP,SS + PP or TT in the SM0M-(7^, y^) scheme. The other SMOM 
schemes also show similarly small mixing coefficients, while the mixing is artificially enhanced 
through the pion pole contribution in the RI-MOM scheme. Since the mixing coefficients are 
found to be at least four orders of magnitude smaller than the multiplicative factor Zj i , we con- 
clude that the mixing can be safely neglected even at the high statistical accuracy reached in our 
computation. In the following we define the renormalization factor for Bk as the multiplicative Z 
factor only. 

3. Step scaling functions 



Following Reference [13811 we can compute the step scaling functions Og^. In this reference a 
comparison of the continuum non-perturbative step scaling functions with the perturbative results 
was proposed as a means to identify the "best" scheme for conversion to MS. It was observed that 
the SMOM(^, ^) scheme agreed very well with the perturbative running. We also find here that 
this scheme has the smallest residual slope in after removing the perturbative running. 



Details of the step scaling scheme can be found in 1)3811 . we briefly summarize them here. Using 
Equation (f89l ) in the chiral limit on each ensemble we have calculated Zg^ (p, a) for p in the range 
2.0 GeV < p < 3.0 GeV. Because of our twisted boundary conditions we have been able to choose 
the same momentum direction consistently. Thus renormalization constants at the same physical 
scale on both lattices have the same Symanzik expansion and we can perform the continuum 
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FIG. 16: Mixing coefficient at j8 = 2.25 for Oi = Ovv+aa and the operators O2 = Ovv-aa, O3 = Oss-pp, 
O4 = Oss+pp and O5 = Ojt- The data shown has been extrapolated to the chiral limit. 



extrapolation of the ratio, 



ZBf.{spo,a) 



where 5 is a scale factor between 1 and 1.5 and po = 2GeV to obtain 

a^Q Zbk[Pq) 



(90) 



(91) 



The present calculation marks an improvement over Reference 113 811 where the determination of the 



lattice spacing was performed using fits to the static potential and was a large source of statistical 
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and systematic error. Here we use the well determined values of the lattice spacing I19I1 on these 
ensembles, which significantly reduces the error. Figure [17] shows the step scaling functions for all 
four SMOM schemes, and we confirm that the SMOM(^, is very well described by perturbation 
theory. This motivates us to use it as our central value. In these plots we use the opposite conven- 
tion to [38] and plot ^^{SGeV) ^^ere s varies between | and 1. The values of aB^(2GeV,3GeV) 
and the corresponding error budgets are presented in Table IXIVl 



scheme 


MOM 


SMOM (7^,7^) 


SMOM (7^,^) 


SMOM (^,7^) 


SMOM (^,^) 


aB^(2GeV,3GeV) 


0.98457 


0.98346 


0.93783 


1.00893 


0.96189 


Stat 


0.00352 


0.00091 


0.00154 


0.00186 


0.00073 


Ms 


0.01041 


0.00075 


0.00382 


0.00056 


0.00012 


V-A 


0.00068 


0.00066 


0.00008 


0.00042 


0.00007 


Total 


0.01101 


0.00135 


0.00412 


0.00199 


0.00075 



TABLE XIV: Scaling factor aBj^(2GeV, 3GeV) from 2 to 3GeV for each scheme. The values are the 
reciprocal of the left most point in Figure [17] The error from the uncertainty in the lattice spacing is now 
folded into the statistical eiTor. 



IV. CHIRAL-CONTINUUM EXTRAPOLATION STRATEGY 

In Reference I19I1 we perform a combined chiral-continuum fit simultaneously to our 1 and 
2 ensemble sets, allowing us to extract the lattice spacing and physical quark masses characteris- 
ing each ensemble set. An ensemble set is a group of ensembles with the same value of /3. When 
extrapolated to physical up/down and strange quark masses, determined via two constraints, we 
determined the lattice spacing of each ensemble set using a third constraint. Thus, with two en- 
semble sets, a total of six constraints are required, and the relation of these constraints between 
the different ensemble sets determines our chosen scaling trajectory to the continuum limit: in 
principle we are free to choose three quantities or ratios as having no corrections in defining our 
scaling trajectory. 

We summarise the chiral-continuum fit procedure and the subsequent determination of the lattice 
scales and physical quark masses below. Throughout we denote masses implicitly shifted by nires 
with a tilde as in m/; these are analogous to a PCAC mass, but as we have good chiral symmetry 
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(e)RI-MOM 

FIG. 17: Continuum limit step scaling functions for all four SMOM schemes (blue) compared with one-loop 
perturbation theory (black). The continuum limit is a simple linear extrapolation in a^. The right, = 1, 
point corresponds to 3 GeV 
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the adjustment is rather small. 



A. Overview of method 

In Reference lioll we simultaneously performed a chiral-continuum fit of the following five quan- 
tities: m^, ruK, mo,, f^^ and fK- After summarising these global fits to obtain lattice spacings 
and quark masses, we will then perform a separate chiral-continuum fit for Bk- We explore two 
alternate sets of fit forms: 

• The first form is obtained through a joint chiral and expansion at next-to-leading order 
in SU{2) chiral perturbation theory (ChPT) and in a^. Throughout our analyses we use 
A;^ = 1 GeV as the chiral scale. For heavy-light quantities such as Bk, rtiK and /k, we 
use SU (2) PQChPT to which the kaon is coupled into the theory at leading order in the 
non-relativistic expansion 

• The second form is obtained from a leading-order analytic expansion about a non-zero un- 



physical pion mass as advocated by Lellouch ^5\, and including corrections. The fit 
forms are linear in the quark masses. By using this approach we lose the ability to take the 
chiral limit and only extrapolate to the non-zero physical point. 



,2 



B. Ideal trajectory to continuum limit 

We must use six quantities to determine the scale, strange mass and the (degenerate) up/down 
mass for each of the two lattice spacings. The discussion can be simplified if we first consider an 
ideal case where we were able to simulate at any quark mass. In this case we would tune the input 
quark masses on both lattices until we obtain m^/ma and mK/ma simultaneously equal to their 
experimentally observed values. 

This would define a non-perturbative, hadronic mass dependent renormalization condition, and the 
freedom we hold in defining the trajectory to the continuum would be absorbed by defining these 
quantities to be artefact free. 
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C. Matching at unphysical quark mass 

In practice, we are not yet able to simulate with the physical quark masses and getting to the phys- 
ical masses involves some degree of interpolation or extrapolation. However, the above strategy 
can be modified to identify the mass parameters for each ensemble which lie on the particular scal- 
ing trajectory by requiring that a pair of mass ratios take on convenient unphysical values rather 
than "real world" observed ratios. 

For example, we can require that the ratios m///m/,/,/, and mhi/nihi^h the values given by one 
pair of input quark masses that were used when generating a particular ensemble. Here the masses 
mil, ifihi and m/,/,/, are the unphysical analogues of m^j:, and m^i for our unphysical choice of m; 
and m/,. Then the pair of matching light and heavy quark masses, (m/, m/,), for a second ensemble 
set with a different value of /3 can be obtained by interpolation in the light quark mass m/. We 
also require a matching value of /n/, on this second ensemble. As we only used one mass value for 
the strange sea quark we apply reweighting to assign the heavy sea quark mass the value This 
self-consistent heavy quark mass reweighting and interpolation to an equal valence mass will be 
performed iteratively. 

We formulate our approach to deal with arbitrarily many /3 values with ensemble set index e. We 
may then define a lattice spacing ratio for each ensemble set e to the primary ensemble set 
1 from the ratio of hhh baryon masses: 



a 



1 



K = = — , (92) 

where this ratio is naturally 1 for e = 1 . 

For the quark masses that yielded matched pseudoscalar and hhh baryon masses we characterize 
the additional logarithmic dependence on a by defining the factors Zf and Zf : 

Zf = (93) 
4 - J^- (94) 



As we approach the continuum limit, standard renormalized perturbation theory implies that phys- 
ically equivalent light and heavy quark masses will be related between two /3 values by the same 
renormalization factor. However, for non-zero lattice spacing we expect Zf ^ Zf . Further as 
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— !■ these factors each approach unity. This implies [jiol that: 

Zl =Zf (l+c„[(aM'-(a')']). (95) 

While the coefficient Cm must vanish as mi ^ m^, we have not written it as proportional to m/, — m/ 
because the low energy matrix elements of the dimension 6 operators which give rise to these 
0{a^) corrections will contain the more complex infra-red quark mass dependence of low energy 
QCD. In fact the difference between these two factors is at or below the 1% level and, as can be 
seen from Table IXVl they were numerically indistinguishable in our study lll9n . Never-the-less we 
treat them as two independent quantities in our fits. 

When performing an extrapolation in quark mass using both of the available ensembles, it is con- 
venient to employ a mass renormalization scheme which is closely related to the mass parameters 
used in those simulations. Thus, for any simulated quark mass on any ensemble set e, we introduce 
an equivalent, matched quark mass , expressed in lattice units on our 1 ensemble set: 

m} = Z} Rl m} for / = Z or h. (96) 

This represents a convenient but unconventional renormalization scheme where Z,„ is defined 
to be unity for our finest lattice spacing. This non-canonical choice of renormalization scheme can 
of course be transformed to MS at a later stage. 

The matching prescription ensures that the trajectory to the continuum is defined such that the 
masses of certain simulated pion-like, kaon-like, and ^2-like particles are lattice artefact free. In 
principle, these states are only lattice artefact free at the specific simulated masses /«/ and m/, used 
to define the fixed factors Z/ and Zf^ in Equation (|96l ). However in some neighbourhood {5,ni-,8,ni^) 
of this simulation point the variations in the factors Z/ and Z/, will be sufficiently small to be 
neglected. Since Z/ and Z/, are already themselves indistinguishable, we can safely neglect the 
variations in Z/ as m/ varies between zero and any of the (0.005, 0.01) and (0.004, 0.006, 0.008) 
quark mass values in our two ensembles. Likewise, we will treat Z/, as constant for 5,,,,^ within 
20% of rrih- Thus, by taking a simulated pion-like object to be artefact free for one of these values 
of mi we can view artefacts in all pions to be small, even in the chiral limit. 

D. SU(2) power-counting 

As in I19I1 we view the light quark mass and expansions as a double power series, and work 
only to NLO in this double series. We choose the quark masses on each ensemble set such that the 
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ratios of some reference pseudoscalar masses to the hhh baryon mass remain fixed. Consider the 
continuum SU(2) expression for the pion mass: 



i^((24^'-4^V2(2L«-Lf [. (97) 

where all quantities are expressed in physical units and 

Xi = 2Bmi (98) 

depends on the definition of the light quark mass m/. When we consider this in an expansion at 
non-zero lattice spacing, we represent B and fhi in our matched lattice scheme as 

2B^ m] 



= TTlf ■ (99) 



As the LEC B is scheme dependent we have used our freedom to define a scheme where it simply 
multiplies the matched bare quark mass on our 1 ensemble. Our matching at non-zero quark 
mass can be introduced to the fit directly with no further a^counter terms as the leading order 
dependence away from our match point has been argued above to be small. For B and m expressed 
in this scheme there are also no order counter terms. 

In fact, we note that if we were to apply Equation (|97l) in independent fits to dimensionless masses 
on each ensemble set, and if the NLO EEC's turned out to be the same (something that our com- 
bined fit constrains to be the case), then our scaling trajectory would require %] to be matched 
in the same way as our earlier matching strategy, that is, xf {(^^ /"^hhh)^ would be required to be 
unchanged along the trajectory. 

These constraints of identical NLO EEC's on both ensembles and fitting our data at the (simu- 
lated) match point would induce the same relation between bare B's on each ensemble that arises 
naturally in our matching approach: 

Xi = {a^ )'^B^ m] = {a"" )'^B' mf (100) 



and thus 



(101) 
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Quantities not used to set quark masses and lattice scales acquire dependence at leading order 
but keep only the continuum portions of next-to-leading order mass-expansion terms. For example, 
the SU(2), partially quenched, light pseudoscalar decay constant for a meson composed of quarks 
with masses m/ and is given by 

At fixed heavy quark mass, we take the partially quenched light quark mass dependence of the 
kaon mass and decay constant as: 

mj, = B"^)!^,,)™,. { 1 + h^x, + M|^);,,| ,103) 



and 



fxh = /(^n'^/Ojl+C^W^'} 



+/(^) (m/0 {+^Xi + ^Xx - ^ 



■log^ + ^^log^^ 



(104) 



These formula have validity once the lattice results have been reweighted so that both valence and 
sea heavy quark masses take the value m/,. 
For the kaon bag parameter we use: 



nxh pO 



1 , 2 , (^OXI . XxC\ X\ 1 / Xx 



E. Analytic expansions 



(105) 



We also consider first order Taylor expansions about a non-zero quark mass m'", in the style of 



14511 ■ By using this approach we lose the ability to take the chiral limit and only extrapolate to 
the non-zero physical point. In fact our ansatz for m;j has a (small when fitted) constant term that 
requires some form of chiral curvature (at smaller masses) to satisfy Goldstone's theorem. Again, 
we apply a power counting rule in a double expansion in 5m and a^. 

For the mass of the pion composed of valence quarks with masses m^, my and as a function of light 
sea quark mass m/ and fixed sea strange mass we write the average valence mass in a meson as 
my = -L^pi. and use the ansatz 

m\ = C'^^ + C'^^ (m, - m'") + C"^'' (m/ - m'") . (106) 
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There is no 0{a^) term at the match point and so no correction to Cq" . Thus within our power 
counting we could equivalently use 

ml = C;,"'' + C'f '^mv + C'^'^m, (107) 

where for convenience we redefine C'^^ between Equations (11061) and (11071) . For decay constants, 
which do not vanish in the chiral limit, the 0{a^) term is not sensitive to the choice of expansion 
point: 

/// = +C/a2] +C('^(mv - m'") +c('^(m/ - m'") (108) 

= C{;'[\+Cfa^]+C{'m^ + C{''m,, (109) 

where again C-[^ has been redefined between Equations (11081) and (11091) . At fixed valence and sea 
strange mass my = m/, = m^, we take the dependence on the light valence quark mass m^ and light 
sea quark mass mi of the kaon mass, kaon decay constant, and kaon bag parameter as 

ml^ = CQ'' + C""'{m^-m"')+C2''{mi-,n"') (110) 
^ C'^'^ + Cl^^m^ + C^'^mi, (111) 

f.h = cl,'[l+Cf,a^]+c{'{m,-m"^)+c('{mi-m"^) 

= €^"[1+ Cf^a^] + Cl^m, + c{^mi , (112) 

Bf = co{l+Caa^)+ci{mi-m"')+Cy{m^-m"') 

= co{l + CaO^) + cifhi + Cyfrix , (113) 

where again the parameters Cq"^, clf and cq have been redefined between each pair of equations, 
and implicitly depend on the strange quark mass. 



V. CHIRAL-CONTINUUM EXTRAPOLATION RESULTS 

In this section we present the joint chiral-continuum extrapolation of our data. 



A. Fitting procedure 

In References J^, we performed correlated fits where the correlation matrix is obtained by 
taking increasing numbers of the leading eigenvectors. We find no significant difference over 
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uncorrelated fit results within our limited ability to estimate the correlation matrix. Hence for this 
analysis and those in References |y,[3] we use uncorrelated fits. 

In order to perform our fits, which include forms valid only for fixed strange mass, we are faced 
with the problem that the physical strange mass is an output of our calculation. Thus the com- 
bined chiral-continuum fit procedure is necessarily iterative. The details of the procedure are 
documented in Reference UM, and it suffices to note here that the iterative process terminates 
when the fixed strange mass forms produce a prediction for that is consistent with the guess 
Ms to which our data was interpolated. When doing this we use reweighting to adjust all pionic 
observables to the current strange mass guess for each ensemble. For kaon and ^ observables a lin- 
ear interpolation between the (unreweighted) unitary measurement, and a second valence strange 
(reweighted-to-be-unitary) measurement suffices to obtain that observable for = m/, = mf 



'guess 



B. Scaling analysis 

As discussed in Section IWl we match our lattice data using ratios of hadronic masses — and 



We choose a specific simulated value of {fhi^rh}^) on the ensemble set M to which the 
other ensemble sets are matched. We refer to this as the match point. The choice of the match 
point defines a particular trajectory along which we approach the continuum limit. Although 
the physical predictions do not depend upon the particular trajectory, certain match points are 
favourable due to the quality of the data at the match point and the range over which the data 
must be interpolated/extrapolated on the other ensemble sets to perform this matching. The ideal 
point has as small a statistical error as possible and lies within the range of simulated data on 
all of the matched ensemble sets such that only a small interpolation is required. In practice, the 
errors on the mass ratios at the match point can be reduced by simultaneously fitting to all partially 

quenched simulated data on the ensemble set M and interpolating to the match point which lies on 

I 

the unitary curve. Further details of the procedure are documented in [ll9]. 

As previously mentioned, the primary ensemble set is chosen to be that with the finest lattice 
spacing; our 32^ x 64, = 2.28 GeV lattice (ensemble 1 ). As we have only one other ensemble 
set, we henceforth drop the superscript on the lattice spacing and quark mass ratios. 



In Table IXVl we give the values [|l9|] for Z/, Z/, and Ra obtained by using several match points on 
both ensemble sets M G {1 ,2 }. Subject to the condition that we require a match point within the 
range of simulated data, we can discard the first and last entries. From the remaining, we choose 
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the values Z, = 0.983(9), Z,, = 0.975(7) and Ra = 0.759(5) from the second entry with M = 1 and 
(m/,m/J^ = (0.006,0.03) as our final values. The consistency is excellent, and these are taken as 
input to our chiral-continuum extrapolation for Bk- 







{miY {mi,Y Zi 


Z/, 


Ra 


A 


0.004 


0.03 


0.00312(13) 0.03804(79) 0.980(15) 0.977(11) 0.7623(71) 


A 


0.006 


0.03 


0.00581(12) 0.03829(51) 0.983(9) 


0.975(7) 


0.7591(46) 


A 


0.008 


0.03 


0.00856(19) 0.03856(63) 0.981(10) 


0.973(8) 


0.7556(58) 


6 


0.005 


0.04 


0.00541(10) 0.03136(48) 0.980(12) 


0.976(8) 


0.7604(55) 


B 


0.01 


0.04 


0.00899(18) 0.03078(56) 0.977(11) 


0.969(9) 


0.7520(69) 



TABLE XV: Values of the quark mass ratios Z/ and Z/, and the lattice spacing ratio Ra determined by match- 
ing at five points over both ensemble sets. Quark masses are quoted without the additive mres correction. 



C. Combined analysis procedure for Bk 



In Reference [|l9h we obtained the the lattice spacings and physical light and strange quark masses 



given in Table IXVTl from our two combined analysis procedures. These are taken as input to our 
fits to Bk in the present calculation. This table also contains the values of the leading-order SU (2) 
ChPT LECs B and / obtainedlH from fitting m^^ and /;;:, and which are used as input to our Bk 
analysis in order to reduce the number of degrees of freedom in the NLO PQChPT fit form. 
In principle, the matrix element fit could be included in our main combined fit analysis, allowing 
these data to constrain the ratio B/ f^. In practice however, this constraint is very weak as com- 
pared to those from iti:^ and fjz, so the Bk analysis can be decoupled from the main analysis. On 
the second line of Table IXVTl we have given the lattice parameters obtained by an NLO PQChPT 
fit with finite volume effects included b y co rrecting the chiral logarithms using the corresponding 



finite volume sum of Bessel functions [|46h . These are propagated through to our analysis of the 
finite volume corrections to Bk- 

Our data are reweighted/interpolated to the physical strange quark mass prior to the fit, as discussed 
above. The data are given in Tables IVl and IVTl We fit this data with both ChPT and analytic forms, 
Equations (11051) and (11131) . fitting the NLO PQChPT form of Equation (11051) both with and without 
finite volume corrections in order to estimate the finite volume systematic error. 
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Fit 








B(GeV) 


/(GeV) 


NLO PQChPT 
NLO PQChPT+FV 
LO Analytic 


2.28(3) 1.73(2) 
2.28(3) 1.73(2) 
2.29(3) 1.74(2) 


0.00099(3) 0.00133(4) 
0.00101(3) 0.00136(4) 
0.00105(6) 0.00140(9) 


0.0278(7) 0.0376(11) 
0.0278(7) 0.0375(11) 
0.0277(7) 0.0374(11) 


4.13(8) 
4.04(7) 


0.107(2) 
0.110(2) 



TABLE XVI: Parameters of the 1 and 2 ensemble sets determined from a combined fit using the fit form 
given in the first column. We also include the LO ChPT LECs B and / that are used to constrain the fits to 
Bk. 

Fit B^(2GeV) Bf^(3GeV) 

NLO PQChPT 0.544(5) 0.523(5) 
NLO PQChPT+FV 0.542(5) 0.521(5) 
LO Analytic 0.557(5) 0.536(5) 

TABLE XVII: B^^(2GeV) as obtained by a combined fit to the data at the physical strange quark mass 
using an NLO PQChPT fit form and a LO analytic fit form. The second line contains the NLO PQChPT fit 
with finite volume corrections included, from which we estimate the finite volume systematic by comparing 
to the fit without coiTcctions. En^ors are statistical only and do not include the error on the renormalisation 
coefficient. 

Note that these equations are applied with strange quark mass fixed to its physical value having 

linearly interpolated and reweighted the data to the physical strange quark mass. 

We renormalize the Bk data using the renormalization constants determined in Section ITlI DI prior 

to performing our fit. Thus the fit is performed seperately for each of the schemes SMOM(^, ^) and 

SM0M(7^, 7^), and for both 2 GeV and 3 GeV matching scales. The central value is taken from 

the SMOM(^,^) scheme, and the SMOM(7^,7^) contributes to determining the renormalisation 

error. 

Performing the fits, we obtain the results given in Table IXVTll where the quoted errors are statisti- 
cal only. Here we have also included an NLO PQChPT fit with finite volume corrections, which is 
used below to estimate the finite volume systematic. The fit parameters are given in Tables IXVlIII 
and IXiXl 

Figure [TS] and [19] display the partially quenched light quark valence and sea mass dependence of 
both our SU(2) and analytic fit forms to kaon matrix element data with one valence quark mass 
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Parameter 


NLO PQChPT 


NLO PQChPT-FFV 




2 GeV 


3 GeV 


2 GeV 


3 GeV 


Ca 
Co 
Cl 


0.533(5) 
0.06(4) 
-0.0060(8) 
0.0061(3) 


0.513(5) 
0.08(4) 
-0.0060(8) 
0.0062(3) 


0.531(5) 
0.05(4) 
-0.0062(8) 
0.0071(4) 


0.511(5) 
0.08(4) 
-0.0062(8) 
0.0071(4) 



TABLE XVIII: Fit parameters of the NLO PQChPT fits to the Bk matrix element, with and without finite 
volume corrections. 



Parameter 


Result 




2 GeV 


3 GeV 


Co 


0.554(5) 


0.534(5) 


Ca 


0.06(4) 


0.08(3) 


Cl 


0.2(3) 


0.2(3) 


Cv 


0.9(1) 


0.9(1) 



TABLE XIX: Fit parameters of the leading order analytic fit to the Bk matrix element. 



set to the physical strange mass, and the sea heavy quark mass reweighted to the physical strange 
mass. Our previous work [j4]] contained small indications in the corresponding plot for curvature 
consistent with NLO ChPT. These have become less pronounced in our doubled data set and also 
not supported by the higher precision data from the second lattice spacing. 

Figure [20] shows the continuum limit chiral extrapolation, overlaid by the data corrected to the 
continuum limit using the fit parameters describing dependence. Figure [21] shows the same fits 
overlaid with the uncorrected data. By comparing these plots, the weak lattice spacing dependence 
of the data is apparent. 



D. Systematic errors on Bk 



Due to our combined analysis technique, and our use of reweighting in the strange sea sector, we 
eliminate systematic errors associated with discretisation effects and the untuned strange quark 
mass that were present in our previous analysis [|3|]. The remaining sources of systematic error are 
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> 

U 

0.565 - 



ra^ = 0.005 
ra =0.01 



0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 0.011 



0.002 0.004 0.006 0.008 0.01 0.012 0.014 



FIG. 18: Partially quenched light valence mass dependence of for the three (32^) 1 ensembles (left 
panel) and two (24^) 2 ensembles (right panel) at a valence strange quark mass fixed to be the physi- 
cal strange mass, and after reweighting in the heavier sea quark mass to the physical strange mass. The 
overlayed curves are the partially quenched SU(2) chiral perturbation theory expressions used in our fits. 




. m^, = 0.004 
. m^, = 0.006 
. m , = 0.008 



0.565 - 




0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 0.011 



0.002 0.004 0.006 0.008 0.01 0.012 0.014 



FIG. 19: Partially quenched light valence mass dependence of Bk for the three (32^) 1 ensembles (left 
panel) and two (24^) 2 ensembles (right panel) at a valence strange quai^k mass fixed to be the physi- 
cal strange mass, and after reweighting in the heavier sea quark mass to the physical strange mass. The 
overlayed hnes represent analytic fits to this data. 



those arising due to the chiral extrapolation, finite volume effects and the renormalization. The 
systematic errors on the renormalization coefficients were discussed in Section Hill We discuss the 
remaining contributions below. 
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0.02 



m, (GeV) 



FIG. 20: The continuum limit chiral extrapolation obtained from our global fits using NLO SU (2) PQChPT 
and LO analytic fits. The data is shown corrected to the continuum limit using the ff{a^) corrections 
obtained from both fit forms. 




♦ 32 data (uncorrected) 

♦ 24^ data (uncorrected) 

— LO analytic 

— NL0SU(2)ChPT 



0.01 



0.02 



0.025 



m, (GeV) 



FIG. 21: The continuum limit chiral extrapolation obtained from our global fits using NLO SU (2) PQChPT 
and LO analytic fits. As opposed to in Figure |20l the data plotted here has not been corrected to the 
continuum limit. The fit curves plotted are those performed to the continuum data as before. 



7. Chiral fit systematics 



n 



In Reference [|4|,ll91] we showed that a continuum fit to our two lattices using NLO SU{2) PQChPT 
fit forms gives a value for f:;i that is ~ 10% too low after finite volume effects are included. 
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Although this is of the magnitude expected for naturally sized NNLO contributions, we show in 
Reference 1I47I1 that a full NNLO fit to our data is heavily dependent on the priors used to constrain 
the fit and thus has little predictive power. We also considered an alternate fit form obtained from 
an analytic expansion at leading order about a non-zero unphysical pion mass, as advocated by 
Lellouch 145]. We are able to fit all of our data successfully, and obtain a result that is much closer 
to the known physical value for /;;:. We observed that the difference between the analytic and the 
ChPT fit results in this case provides a good estimate of the systematic error associated with the 
chiral fit form] 18, 3]- We concluded that comparing ChPT and LO analytic fits is likely a good, 
robust method of estimating the systematic error for other quantities such as Bk- Both approaches 
must converge upon the physical value as the simulated quark masses approach the physical point. 
The result of the LO analytic fit to Bk is given alongside the NLO PQChPT results and those 
with NLO PQChPT including finite volume effects in Table IXVIII To combine these in a final 
prediction, we follow jiol and note that both the analytic and finite volume NLO PQChPT fits 
are reasonable extrapolation methods that can be justified in distinct limiting cases: the analytic 
form is certainly the correct approach when we have data sufficiently close to the physical point 
regardless of whether we are in the chiral regime, while the NLO form including finite volume 
effects is also certainly correct when the data and physical point lie within the chiral regime. 
Given our experience with fji, and following the approach taken in Ill9ll we take our central value as 
the average of those obtained with the analytic extrapolation form, and the finite volume corrected 
SU(2) NLO forms. We take the difference between these to estimate a chiral fit systematic error as 
(AB/f);^ = 0.014 (2.6%). We take the full difference as the systematic and believe this is a prudent 
and conservative approach. 

Another reasonable data driven method would take half the difference as the error estimate; this 
would assume that the analytic extrapolation is a hard upper bound on the mass dependence, and 
that the NLO form is a hard lower bound - given the flexibility in unconstrained NNLO ChPT 
forms this would appear to be too optimistic. 

We also note that within the mass range of the data our SU(2) NLO fit estimates the biggest 
correction to be around 8% of the value in the two flavor chiral limit (0.56 vs 0.517). Squaring 
this term would suggest a naive estimate of NNLO effects at around 0.5%, which is substantially 
below our more conservative chiral extrapolation error. 
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2. Finite volume systematics 

We estimate finite volume corrections to our result from finite-volume PQChPT. As shown in 
Reference 141] these corrections are obtained from the standard PQChPT forms by replacing the 
NLO chiral logarithms with sums over modified Bessel functions of the second kind. 
The result for this fit is given in Table IXVIIl Comparing this to the uncorrected result we estimate 
a finite volume error of {1^k)vn = 0.002 (0.4%). 

E. Continuum prediction for Bk 

Combining our central value and the systematic uncertainties discussed above, we quote a predic- 
tion for Bfc using either the = jl^ = (2 GeV)^ renormalization scale, 

5^S(2GeV)=0.549(5)stat(15);^(2)Fv(21) NPR- (114) 
or the = jl^ = (3 GeV)^ renormalization scale 

5fS(3GeV) =0.529(5)stat(15);^(2)Fv(ll)NPR. (115) 

The latter is our preferred central value as our systematic error for the renormalization is halved. 
This can be converted to the common RGI scheme for comparison and phenomenological appli- 
cation: 

' ' (116) 



5|Gi = o.749(7)stat(21);t(3)Fv(15)NPR, 



and adding all sources of error in quadrature we obtain 

' (117) 



5|GI = 0.749(27)combined, 



corresponding to an overall error of 3.6%. 



VI. CONCLUSIONS 



In this paper we have calculated Bk to 3.6% precision with 2-1-1 flavours of dynamical quarks and, 
for the first time, in the continuum limit with a lattice action with good chiral symmetry. The result 
is presented in Equation (|116l) (or equivalently in (11171) ). 

Our calculation of this important quantity has exploited several significant improvements in lattice 
techniques which we have been developing for more than a decade. These include: a) the use of 
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Publication 



/ 



B 



RGI 



This work 



2+1 0.749(7)(26) 



Bae'10[52] 2+1 0.724(12)(43) 

RBC-UKQCD'09[18J 2+1 0.737(26) 

Aubin'09 [53J 2+1 0.724(8)(29) 

RBC-UKQCD'07L3] 2+1 0.720(13)(37) 

ETMC'10[54J 2 0.729(30) 

ETMC'09 [55J 2 0.73(3)(3) 

JLQCD'08 [56J 2 0.758(6)(71) 

TABLE XX: A comparison of our result for Bk with those of other recent calculations with dynamical 
fermions. Here / denotes the number of dynamical quark flavours. Where separate errors are quoted, the 
first enw is statistical and the second is systematic. 



domain wall fermions with good chiral symmetry [[g, |48D, b) the implementation of domain wall 



TM2A 



49 -E 



^914510, and c) 



fermions in dynamical simulations with 2 + 1 flavours of light quarks |13|, 
the use of SU(2) ChPT for chiral extrapolations of 2+1 flavour simulations, first exploited by the 
RBC-UKQCD coUaborations UM- 

The present calculation of Bj^ includes a particularly careful treatment of the renormalization. We 
have introduced several new momentum renormalization schemes (based on the original works of 
I26I1 and of [3^ as explained in detail in Section lllll). and our renormalization also includes, for 
the first time, the improved scaling procedure of 



38h. 



The small increase in our central value for Bk in this work and in jlSll compared to yl, IJ] has 
arisen partly from significant improvements in our approach to renormalization as well as from 
taking the continuum limit. The difference is within the previously budgeted errors for these 
sources, and a large component of this small shift arises from taking the central value from a new, 
non-exceptional momentum scheme using the perturbative results derived in this paper. 
Our result for Bk is compared to other recent calculations in Table IXXl Since all the results in this 
table, except for those of Reference [|52|] and the current work, used the original RI-MOM scheme, 
there is a substantial correlation in the perturbative systematics between these five calculations. 
Thus the additional renormalization schemes introduced in this paper give added confidence to the 
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estimates of the systematic error from this source. 

In the remainder of this section we briefly discuss the significance of the recent lattice results for 
Bk and the prospects for improving the precision still further. 



A. Significance of lattice results of Bk 

Flavour physics will continue to be central to the exploration of the limits of the standard model, to 
searches for new physics and to the eventual understanding of the fundamental theoretical frame- 
work of physics beyond the standard model. An important tool in this endeavour is the interpre- 
tation of experimental data in terms of the unitarity triangle where, in general, the remarkable 
consistency of the information from different processes places significant constraints on the pos- 
sible parameter space of new models. Having said this^number of tensions have arisen in recent 



years; possible inconsistencies at a 1.5 — 3c7 level [|57l- l60(l which certainly merit further investiga- 
tion. The lattice results for Bk contribute to these tensions as we now briefly explain. 
Lattice calculations are necessary to evaluate the hadronic effects in tests of the unitarity of the 
CKM matrix and our results for Bk, used in conjunction with the experimental determination of 
£k, the indirect CP violation parameter monitoring — nn, are a major ingredient in tests of 
the CKM paradigm (see Equation Qi). We illustrate this here with one example, exploiting lattice 
inputs not only for Bk but also for the semileptonic B ^ K,p and B D,D* formfactors (used 
to determine Vub/Vcb) and the SU(3) breaking ratio, ^, which contains the hadronic effects in 
the ratio of the mixings of B^ mesons and B^ mesons. With these three key lattice inputs a nice 



prediction, sin 2/3 = 0.75 ±0.04 11571 - 15 911 . emerges. This can be compared with direct experimental 
measurements from the time-dependent CP asymmetry in the golden mode, Bj J /W^s which 
gives, sin 2/3'^/'^^' = 0.681 ±0.025 iJ, which is within 2a of the Standard Model prediction with 



6li 



6211 who stress the need to 



the lattice input. A similar tension is found in References [|5, 
include better approximations to the theoretical expression for Ck now that Bk is known to such 
good precision. These improvements include terms proportional to ImAo/ReAo (where Aq is the 
K — )• KK amplitude with the two pions in a state with isospin 0) and the recognition that the phase 
arctan(2AM/f/Ar) is not precisely equal to ;r/4 {/SMk and AT are the differences of the masses 
and widths of the Kl and ^5 mesons). 

From the above discussion it is clear that lattice calculations of weak matrix elements in general, 
and of Bk in particular, in conjunction with experiments, are providing ever more precise tests 
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of the CKM explanation for CP violation. Of course our ambitions do not stop here; even if 
the small tension between the Standard Model prediction for sin(2/3) and its direct determination 
disappears on closer scrutiny, the 0(10%) difference in the central values still leaves ample room 
for new physics which we wish to squeeze still further. In the next subsection we discuss the 
prospects for improved precision in the determination of Bk and of course it must be remembered 
that improvements in the determination of other inputs, including and Vcb will also be necessary 
(recently it was shown that the use of V^^ with its significant error, can be replaced by information 
from the leptonic 5 — >■ TV branching ratio and lattice results on the decay constant fgj and the 



mixing parameter 1631] ) 



B. Prospects for Bk with one percent scale precision 

It is interesting to analyse our error budget and to assess what future gains in precision can be 
made in the determination of Bk- In particular, we consider here what would be required to obtain 
Bk with one percent scale precision. 

Currently, our dominant uncertainty is the 3% error arising from the chiral extrapolation. This will 
be addressed by simulations at or near the physical quark masses, some of which are presently be- 
ing undertaken by RBC and UKQCD. Although expensive, these are affordable, even with current 
computer technology. We can therefore envisage these to be under control at the one percent level 
in a few years. 

The 2% renormalization error is partly associated with the low scale at which we presently apply 
one-loop matching and two-loop running to our operators. This uncertainty can be reduced in two 



ways: firstly the scale can be raised at modest expense using a step scaling technique 13 811 . perhaps 
raising the matching scale from around 3 GeV to approximately 10 GeV, reducing the error on 
our one-loop matching from 2% to around 1%. A larger gain would be obtained by extending 
the perturbative calculations presented in this paper to the next order, leading to an expected a| 
error of around 0.7%. The gain from step scaling is of course increased by higher order match- 
ing, and one might expect a step scaled matching to attain 0.2% renormalization precision for an 
renormalization error. Such a two-loop calculation has been performed for the determination 
of light-quark masses jsS, 34] contributing to the improved lattice determination of these quanti- 



ties uM- Given the importance of a precise determination of Bk, we would hope and expect that 



the two-loop matching calculation will be performed soon. 
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The remaining statistical and finite volume errors are small, and not unduly expensive to reduce 
still further as this increases computational cost by only modest factors. 

We conclude therefore that we can expect to determine Bk at the one percent scale over the next 
few years. What is perhaps more challenging is for lattice simulations to contribute in other ways 
to the determination of subdominant corrections to the theoretical expression for £k, for example 
the long-distance contributions and the direct computation of ^ — )■ nn decay amplitudes; the status 



of our endeavours in this direction are summarised in [ 16l . 
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